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Preface 


Panel  flutter,  an  aeroelastic  stability  structural  problem,  has  been  a  research  topic  for  more  than 
three  decades.  When  a  flight  vehicle  flies  at  a  supersonic  speed  in  the  air,  some  skin  panels  may 
experience  high  level  vibrations  and  fail  due  to  the  aerodynamic  pressure  on  the  vehicle  surface.  This 
aeroelastically  induced,  self-excited  motion  has  been  described  as  panel  flutter.  In  the  early  stage  of 
this  research,  before  1960s,  it  was  usually  assumed  that  the  phenomena  could  be  adequately 
understood  in  terms  of  linear  models  and  inviscid  aerodyanmics  and  that  the  questions  of  self-excited 
stability  and  response  to  external  disturbances  could  be  readily  separated.  Later,  from  the  late  of  1960s 
and  early  1970s,  it  has  been  found  that  these  assumptions  are  invalidated  most  frequently  for  the 
aeroelastic  behavior  of  plates  and  shells.  Great  progress  has  been  made  in  understanding  the  effects  of 
nonlinearities,  viscous  flow,  and  complex  relation  between  the  stability  and  response  since  then.  Most 
recently,  the  smart  material  and  structural  technique  has  been  implemented  into  this  area  to  suppress 
the  flutter  of  panels. 

However,  most  of  the  researches  on  this  area  are  concentrated  on  the  stmctural  side,  i.e.,  panel  or 
plate.  In  these  researches,  the  approximate  theories,  such  as  qusi-steady  piston  theory,  full  linearized 
(inviscid)  potential  flow  theory,  etc.,  are  used  to  estimate  the  aerodynamic  pressure.  This  kind  of  linear 
aerodynamics  may  not  be  adequate  to  predict  the  dynamic  characteristics  of  the  fluid  and  structure 
because  the  fluid  flow  is  strongly  nonlinear  at  the  transonic  and  supersonic  speeds.  As  we  know,  the 
high-fidelity  equations,  such  as  Euler  or  Navier-Stokes  equations,  can  predict  the  flow  characteristics 
more  accurately.  One  of  the  important  reasons  that  the  high-fidehty  equations  have  not  been  used  to 
predict  the  aerodynamic  loads  is  that  the  corresponding  niunerical  simulation  is  very  computational 
expensive.  With  the  fast  development  of  the  computer  techniques,  the  full  analysis  of  the  nonlinear 
panel  flutter  coupled  with  the  Euler  or  Navier-Stokes  flow  equations  becomes  possible. 

The  Air  Force  Office  of  Scientific  Research  (AFOSR),  Department  of  Defense  (DOD),  sponsored 
a  grant  (F49620-98-1-0396)  through  the  University  of  Arkansas  at  Fayetteville  during  the  period  of 
April  1st  1998  and  March  31st  2001.  The  development  of  an  aeroelastic  solver  which  combined  the 
Computational  Structural  Dynamics  (CSD)  and  the  Computational  Fluid  Dynamic  (CFD)  was  required 
in  this  grant.  Such  a  code  is  able  to  predict  the  change  in  shape  or  position  of  a  stmcture  due  to  a 
calculated  fluid  pressine  exerted  on  it,  and  able  to  model  the  flow  of  fluid  around  the  stmcture.  As 
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required,  the  aeroelastic  solver  has  been  developed.  Some  important  and  usefiil  results  have  been 
proposed. 

This  report  sets  out  to  explain  the  principles  used  in  the  aeroelastic  solver.  In  the  structural  side, 
two  formulations,  total  Lagrangian  and  co-rotational  formulations,  for  the  nonlinear  analysis  of  the 
thin  plate  are  discussed  in  detail.  In  the  fluid  side,  the  finite  volume  method  has  been  used  to  solve  the 
Euler  and  Navier-Stokes  equations.  To  accurately  model  the  problems  with  moving  and  deforming 
boundaries,  the  deforming  mesh  schemes  are  provided  and  discussed  to  make  the  fluid  grid  conform  to 
the  changing  shape  of  the  boundary.  The  numerical  simulation  of  the  aeroelastic  solver  is  provided 
using  a  specific  panel.  Several  issues  of  the  panel  flutter  are  surveyed  at  the  end  of  this  report. 

We  would  like  to  express  our  special  thanks  to  the  AFOSR,  DOD,  for  their  sponsorship  of  this 
research.  Thanks  also  go  out  to  Dr.  R.  E.  Gordnier  and  Dr.  M.  R.  Visbal  at  the  Air  Force  Research 
Laboratory  for  their  help  during  the  research.  We  also  would  like  to  record  our  appreciation  of  many 
helpful  information  and  comments  made  by  the  graduate  students,  Jerry  Scott  R  Prescott,  Yangki 
Jung,  Douglas  Benton,  in  this  laboratory.  Finally,  our  wives  and  children  are  certainly  not  to  be 
omitted  from  our  acknowledgement. 


R.  Parmeer  Selvam  Zu-Qing  Qu 

Fayetteville,  Arkansas,  USA 
November  2001 
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0.1  Importance  of  Current  Panel  Flutter  Analysis 

With  the  resurgent  interest  in  flight  vehicles  such  as  the  High-Speed  Civil  Transport  (HSCT),  the 
X-33  Advanced  Technology  Demonstrator,  the  Reusable  Launch  Vehicle  (RLV),  the  Joint  Strike 
Fighter  (JSF)  and  the  X-38  Spacecraft  using  a  lifting-body  concept  that  will  operate  at 
supersonic/hypersonic  Mach  numbers,  the  need  for  panel  flutter  analysis  has  received  broad 
acknowledgement' . 

The  linear  and  nonlinear  analysis  of  the  panel  flutter  has  been  studied  extensive  dining  the  past 
two  decades'  *.  However,  most  of  the  researches  on  this  area  are  concentrated  on  the  structural  side, 
i.e.,  panel  or  plate.  In  these  researches,  the  approximate  theories,  such  as  qusi-steady  piston  theory’’*, 
full  linearized  (inviscid)  potential  flow  theor/’'",  etc.,  are  used  to  estimate  the  aerodynamic  pressure. 
This  kind  of  linear  aerodynamics  may  not  be  adequate  to  predict  the  dynamic  characteristics  of  the 
fluid  and  stmcture  because  the  fluid  flow  is  strongly  nonlinear  at  the  transonic  and  supersonic  speeds. 
As  we  know,  the  high-fidelity  equations,  such  as  Euler  or  Navier-Stokes  equations,  can  predict  the 
flow  characteristics  more  accurately.  One  of  the  important  reasons  that  the  high-fidelity  equations 
have  not  been  used  to  predict  the  aerodynamic  loads  is  that  the  corresponding  numerical  simulation  is 
very  computationally  expensive.  With  the  fast  development  of  the  computer  techniques,  the  full 
analysis  of  the  nonlinear  panel  flutter  coupled  with  the  Euler  or  Navier-Stokes  flow  equations  becomes 
possible. 
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0.2  Description  of  Current  Panei  Flutter  Anaiysis 

Panel  flutter  in  the  supersonic  flow  falls  in  the  category  of  self-excited  or  self-sustained 
oscillations.  Historically,  flutter  and  other  problems  involving  a  fluid  and  a  structure  were  predicted 
using  wind  tunnels.  This  type  of  testing  required  building  small  test  models,  which  are  expensive  and 
time  consuming.  In  addition  to  physical  models  being  costly,  testing  facilities  such  as  wind  tuiuiels 
capable  of  producing  air  speed  beyond  the  speed  of  sound  are  very  rare.  However,  within  the  past 
couple  of  decades,  advances  in  computer  technology  has  greatly  lowered  the  time  and  cost  of  design 
by  enabling  us  to  model  complex  engineering  problems  on  computer  rather  than  building  costly 
physical  test  models.  Through  computer  modeling  it  is  possible  for  engineers  to  predict  the  motion  of 
complex  shapes  in  a  fluid  field  and  correct  failures  and  other  problems  of  a  design  before  the  first 
prototype  is  built. 

The  aerospace  industry  was  among  the  first  to  use  and  further  develop  computer  modeling  to 
predict  all  types  of  phenomenon,  several  of  these  are  dynamic  loads,  moments,  fluid  pressures,  shock 
waves,  and  turbulence.  Initially,  computer  modeling  of  aircraft  was  separated  into  two  fields. 
Computational  Structural  Dynamics  (CSD)  and  Computational  Fluid  Dynamics  (CFD).  In  CSD  a 
structure  is  analyzed  by  assuming  a  fluid  pressiue  for  a  given  air  speed  as  the  structural  load.  Then 
apply  it  to  the  structure  to  obtain  the  displacements  and  stress.  For  CFD  the  flow  of  fluid  is  analyzed 
arormd  an  immobile,  infinitely  stiff  structure. 

As  technology  progressed  to  produce  faster  aircraft  so  did  the  need  to  imderstand  the  nature  of 
fluid  and  its  interaction  with  a  structure  and  more  importantly  to  predict  it.  Highly  maneuverable 
aircraft  such  as  fighters  can  experience  very  complex  fluid  flow  such  as  vortices,  flow  separation, 
shock  waves,  flutter,  and  other  high  speed,  nonlinear  flow  phenomenon. 

Such  problems  create  the  need  to  combine  CFD  and  CSD.  The  unity  of  these  two  programs  is 
known  as  an  aeroelastic  code.  Such  a  code  is  able  to  predict  the  change  in  shape  or  position  of  a 
structure  due  to  a  calculated  fluid  pressure  exerted  on  it,  and  able  to  model  the  flow  of  fluid  aroimd  the 
structure.  However,  neither  CFD  nor  CSD  has  been  perfected  to  give  exact,  computationally  efficient 
results.  Any  errors  produced  in  the  structural  calculations  are  carried  over  to  the  fluid  calculations  and 
vice  versa  via  the  fluid  structure  interaction. 

Fluid  structure  interaction  is  the  exchange  of  information  or  data  from  one  medium  to  the  other. 
This  takes  place  where  the  boundary  of  the  fluid  meets  the  boundary  of  the  structure.  The 
computational  result  of  one  medium  becomes  the  initial  boundary  conditions  of  the  other.  In  a  physical 
problem  the  effect  of  one  particle  of  fluid  on  a  structure  will  change  the  motion  of  the  stracture  and  the 
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fluid  will  react  to  that  change.  In  nature  this  process  happens  almost  instantaneously.  To  accurately 
model  this  behavior  the  equations  of  motion  for  the  fluid,  structure,  and  the  fluid-structure  interaction 
have  to  be  solved  simultaneously.  An  integrate  algorithm  in  which  the  fluid  and  structure  are  modeled 
as  a  single  dynamic  system'*  may  satisfy  this  requirement.  This,  however,  requires  the  development  of 
an  entirely  new  solver.  One  commonly  used  scheme  is  that  the  effects  of  one  medium  on  the  other  is 
delayed  or  lagged  in  time  and  the  results  of  one  medium  can  be  used  to  calculate  the  results  of  the 
other.  The  amount  of  lag  or  time  step  size  affects  the  accuracy  of  the  results.  To  reduce  this  error  the 
Newton-type  subiteration'^  may  be  used. 

Another  obstacle  to  be  overcome  in  the  fluid-structure  interaction  is  the  coordinate  system  that 
describes  the  motion  of  a  fluid  or  a  structure  (The  computational  grid  of  the  flxxid  or  the  structure 
represents  the  coordinate  system).  There  are  two  types  of  coordinate  systems,  the  Euler  and  the 
T  jigrangian  Generally,  the  motion  of  a  fluid  is  described  by  the  Euler  coordinate  system,  where  the 
coordinate  system  does  not  move  with  the  material  particles  in  time.  The  motion  of  a  stracture  is 
described  by  the  T.?igrangian  coordinate  system,  where  the  coordinate  system  moves  with  the  material 
particles  in  time.  When  dealing  with  a  stmcture  that  is  moving  within  a  fluid,  the  fluid  grid  must  be 
updated  at  each  time  step  to  conform  to  the  new  position  of  the  structure.  To  solve  this  issue  two 
things  must  be  taken  into  account,  how  to  modify  the  fluid  equations  to  allow  for  movement  and  still 
get  accurate  results  and,  how  to  create  a  gnd  to  conform  to  the  new  shape  of  the  structure. 

Small  amplitude  linear  structural  theory  indicates  that  there  is  a  critical  dynamic  pressure  above 
which  the  panel  motion  becomes  unstable  and  grows  exponentially  with  time,  hi  reality,  however, 
geometrically  nonlinear  effect  is  present  due  to  vibration  with  large  amplitudes.  The  panel  not  only 
bends  but  also  stretches  with  in-plane  tensile  stresses.  Such  in-plane  tensile  stresses  provide  a 
stabilizing  effect  that  restrains  the  panel  motion  to  bound  limit  cycle  oscillations  with  increasing 
amplitude  as  the  dynamic  pressure  increases.  The  skin  panels  of  the  flight  vehicle  can  thus  withstand 
velocities  beyond  those  predicted  based  on  the  small  amplitude  linear  structural  theory.  For  a  more 
thorough  imderstanding  and  more  realistic  assessment  of  the  panel  flutter  problems,  the  geometrically 
nonlinear  effect  due  to  the  large  amplitude  oscillations  should  be  considered. 


0.3  Organization  of  this  Report 

Geometric  nonlinearity  is  increasingly  being  considered  in  structural  analysis,  especially  for  those 
structures  that  undergo  large  displacement  with  very  little  actual  staining  occurring.  There  are  many 
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woiics  on  the  way  in  which  the  geometrically  non-linear  problems  may  be  formulated.  Three  main 
methods  have  been  currently  employed,  namely,  the  Total  Lagrangian  (TL)  formulation,  the  Updated 
Lagiangian  (UL)  formulation  and  the  Co-Rotational  (CR)  formulation.  Each  of  these  formulation  has 
its  inherent  drawbacks,  merits  and  advocates.  The  details  of  the  pros  and  cons  of  the  formulations  were 
surveyed  using  a  cantilever  beam  in  our  laboratory‘s.  Due  to  the  slow  convergence  of  the  UL 
formulation,  the  TL  and  CR  formulations  will  be  discussed  in  this  report. 

The  large  deflection  of  thin  plates,  which  considers  the  coupling  effect  between  the  flexure  and  the 
in-plane  deformation,  is  investigated  in  Chapter  1  by  using  the  fimte  element  method.  The  TL 
formulation  in  which  the  original  imdeformed  configuration  is  taken  as  the  reference  configuration  is 
to  be  used.  Based  on  the  von-Karman  nonlinear  strain-displacement  relation,  the  dynamic  equations  of 
large  deflection  plates  are  derived  at  first.  The  formulation  of  Discrete  Kirchhoff  Theory  (DKT)  for 
triangular  elements  is  then  presented  for  the  discretization  of  the  thin  plate.  An  iterative  method 
between  the  equations  of  the  flexure  and  membrane  is  used  to  solve  the  nonlinear  equations  at  each 
time  step.  To  make  the  method  efficient,  a  relaxation  factor  was  also  used.  TheNewmark-  method  is 
used  to  approximate  in  time.  Several  numerical  examples  are  also  included  in  this  chapter  to 
demonstrate  the  efficiency  of  the  method  for  nonlinear  dynamic  problems. 

The  CR  formulation  makes  use  of  both  the  original  configuration  and  the  current  deformed 
configuration  to  formulate  the  system  matrices.  Although  restricted  to  small  rotation  between 
iterations,  it  exhibits  good  rate  of  convergence  and,  thereby,  reduces  the  computation  time.  In  Chapter 
2,  CR  formulation  for  static  and  dynamic  analysis  of  plate  structures  using  DKT  triangular  elements  is 
presented.  A  number  of  numerical  examples  are  presented  to  fully  test  the  capabilities  of  this 
formulation.  The  CR  formulation  code  developed  in  this  study  can  be  used  for  more  understanding  and 
realistic  assessment  of  the  panel  flutter  problems. 

As  mentioned  above,  the  fluid  equations  are  usually  derived  using  the  Eulerian  motion  description 
and  the  structural  equations  derived  using  the  Lagrangian  description.  To  accurately  model  problems 
with  moving  and  deforming  boundaries  requires  that  the  fluid  grid  conform  to  the  changing  shape  of 
the  boundary.  To  accommodate  the  movement  of  the  fluid  grid  either  a  new  grid  must  be  generated  at 
each  time  step  (dynamic  regridding)  or  the  differential  equations  for  the  fluid  must  be  modified  to 
allow  the  movement  of  the  existing  fluid  grid  (ALE,  Corotational,  Space-Time).  Once  a  method  is 
found  for  modifying  the  fluid  equations  to  allow  for  grid  movement  a  method  must  be  chosen  to  move 
or  update  the  fluid  grid  to  the  new  position  of  the  stmcture.  In  Chapter  3,  different  methods  used  to 
modify  the  fluid  equations  to  allow  the  movement  of  a  structure  within  the  fluid  domain  are  described. 
For  the  ALE  method,  various  methods  used  to  update  or  move  the  fluid  grid  to  conform  to  a  changing 
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structural  borjndary  are  also  described  and  evaluated.  The  details  of  the  implementation  of  the  Trans- 
Finite  Interpolation  scheme  for  the  deforming  mesh  is  provided  in  this  chapter.  Numerical  simulation 
is  also  included  to  show  the  features  of  this  scheme. 

In  the  fluid  side  of  the  fluid-structure  interaction,  the  fluid  dynamic  computation  must  have  the 
fidelity  to  capture  the  relevant  flow  features  and  provide  accurate  aerodynamic  loads  on  the  structure 
in  developing  a  perfect  fluid-structure  interaction  solver  in  computational  aeroelasticity.  A  finite 
volume  method  is  used  in  Chapter  4  to  establish  computational  solvers  for  Euler  equation  and  Navier- 
Stokes  equation. 

A  numerical  simulation  of  the  three-dimensional  nonlinear  panel  flutter  at  supersonic  flow  is 
performed  using  the  high-fidelity  flow  equations  -  Euler  equations  in  Chapter  5.  These  flow  equations 
are  solved  using  the  Beam-Warming,  alternate-direction,  implicit  scheme.  The  governing  equations  of 
the  thin  panel  are  based  on  the  von-Karman  large  deflection  plate  theory.  Since  the  finite  difference 
method  is  much  more  efficient  than  the  finite  element  method  for  a  simple  structure,  the  former  is  used 
to  discretize  these  governing  equations  in  space  and  Newmark-  P  integration  scheme  is  applied  to 
solve  them  in  time  domain.  The  Newton-like  subiteration  is  implemented  to  eliminate  the  lagging 
errors  associated  with  the  exchange  of  the  pressure  and  deformations  between  the  fluid  and  panel  at 
their  interface.  A  simply  supported  square  panel  is  considered  at  this  time.  The  results  at  the 
supersonic  flow  are  compared  with  those  fi-om  the  aerodynamic  approximations,  such  as  potential  flow 
theory  and  qusi-steady  supersonic  theory.  The  flutter  frequency  and  the  effects  of  the  panel  thickness 
are  also  discussed  in  this  chapter. 
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CHAPTER  1 


NONLINEAR  PLATE  ANALYSIS  USING  TOTAL 
LAGRANGIAN  FORMULATIONS 


1.1  Introduction 

1.2  Finite  Element  Formulation  of  Large  Deflection  Plates 

1.3  DKT  Triangular  Element 

1.4  Static  Nonlinear  Problems  of  Large  Deflection 

1.5  Dynamic  Nonlinear  Problems  of  Large  Deflection 

1.6  Non-Dimensional  Form  of  FEM  Formulations 

1.7  Summary 


1.1  Introduction 

Panel  flutter,  an  aeroelastic  structural  problem,  has  been  a  research  topic  for  the  past  three  decades 
and  has  recently  received  great  interest  When  a  vehicle  flies  at  a  supersonic  speed  in  the  air,  some 
skin  panels  may  experience  high  level  vibrations  and  fail  due  to  the  aerodynamic  pressure  on  the 
vehicle  surface.  This  aeroelastically  induced,  self-excited  motion  has  been  described  as  panel  flutter. 
Experiments  showed  that  there  are  critical  dynamic  pressures  in  panel  flutter.  Below  these  critical 
pressures  the  panel  has  a  random  oscillation  with  small  amplitude  which  is  a  small  fraction  of  the 
panel  thickness.  The  predominant  frequency  components  are  observed  to  be  near  the  lower  panel 
natural  frequencies.  Basically,  the  panel  is  undergoing  a  linear  oscillation.  These  critical  dynamic 
pressures  are  also  called  the  flutter  boundary.  Beyond  this  boundary,  the  amplitude  of  the  panel 
oscillation  grows  rapidly  to  the  order  of  the  panel  thickness.  For  this  case,  the  responses  predicted 
using  the  small  deflection  linear  theory  are  no  longer  applicable,  and,  therefore,  nonlinear  theory 
taking  account  the  effects  of  large  amplitude  should  be  applied. 

The  analytical  solution  of  the  vibrations  of  large  deflection  thin  plates  was  surveyed  by  many 
researchers  more  than  20  years  ago.  Most  of  the  studies,  however,  have  been  concerned  with  circular 
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or  rectangular  plates  due  to  the  difficulty  of  the  mathematical  treatment.  With  the  increased  use  of  the 
plate  in  many  optimum-  or  minimum-weight  designed  built-up  structures,  many  numerical  methods, 
such  as  finite  difference  method*,  the  finite  element  methocf  *,  finite  strips  methocf  ■‘®,  finite  element  - 
transfer  matrix  method"  and  so  on,  were  used  to  discretize  the  plates  in  physical  space. 

The  nine-parameter  bending  plate  element  is  a  simple  element  that  is  widely  used  in  engineering. 
Many  such  triangular  elements  have  been  developed.  They  include  the  element  BCIZ",  discrete 
Kirchhoff  element  DKT",  natural  formulation  element  TRUNC",  free  formulation  element  Tja", 
hybrid  stress  element  etc.  Recently,  a  set  of  three  good  triangular  elements'*  and  a  refined 

direct  stiffness  method  (RDSM)”  have  been  proposed. 

The  formulation  of  elements  based  on  the  discrete  Kirchhoff  theory  for  thin  bending  plates  is 
obtained  by  considering  a  theory  of  plates  including  transverse  shear  deformation.  In  this  case,  the 
independent  quantities  are  the  deflection  and  rotations  and  only  Cf*  continuity  requirements  need  to  be 
satisfied.  The  transverse  shear  energy  is  neglected  altogether  and  the  Kirchhoff  hypothesis  introduced 
in  a  discrete  way  along  the  edges  of  the  element  to  relate  the  relations  to  the  transverse  displacements. 

Batoz^®’^'  studied  the  DKT  element  and  concluded  that  it  is  one  of  the  efficient,  cost  effective  and 
reliable  elements  of  its  class  for  static  bending.  The  convergence  properties  of  the  DKT  element  do  not 
deteriorate  with  an  increase  in  the  element  aspect  ratio  and  which  is  not  so  far  for  other  elements. 
Ratoz^'  also  presented  the  stiffness  matrix  of  the  DKT  element  in  explicitly  form  in  local  coordinate 
system.  This  means  that  the  stiffness  matrix  is  to  be  transformed  to  a  global  coordinate  system  before 
assembly.  A  simple  algorithm  was  proposed  by  Jeyachandrabose  and  Kirkhope“  to  directly  obtain  the 
stiffness  matrix  in  the  global  coordinates.  The  relationship  of  a  series  of  elements  based  on  Reissner- 
Mindlin  assumptions  and  using  discrete  (collocation  type)  constraints  was  discussed  by  Zienkiewicz 
and  Taylor  et  aP.  Most  recently  many  researchers^'’’^  have  been  trying  to  improve  the  accuracy  and 
efficiency  of  the  DKT  triangular  and  rectangular  elements. 

In  this  report,  the  large  deflection  of  thin  plates,  which  considers  the  coupling  effect  between  the 
flexure  and  the  in-plane  deformation,  is  investigated  by  using  finite  element  method.  Based  on  the 
von-Karman  nonlinear  strain-displacement  relation  and  the  DKT  of  triangular  elements,  the  dynamic 
equations  of  large  deflection  plates  are  derived  using  the  finite  element  method.  The  Newmark-  ^ 
method  is  used  to  solve  the  dynamic  equations  of  equilibrium  in  time  domain.  Several  numerical 
examples  are  also  included  to  demonstrate  the  efficiency  of  the  method  for  nonlinear  dynamic 
problems. 
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1.2  Finite  Element  Formulation  of  Large  Deflection  Plates 

1.2.1  Assumptions 

The  following  assumptions  will  be  used  during  the  derivation  of  the  finite  element  formulation  of 
the  large  deflection  plates: 

•  Isotropic  material  obeys  Hook’s  Law  (small  strain); 

•  Thin  plate  (a/t>20),  and  <7^  =  0 ; 

•  Inplane  inertia,  rotatory  inertia  and  transverse  shear  deformation  effects  are  negligible; 

•  Von  Karman  large  deflection  strain-displacement  relations  are  valid.  The  displacements  are  not 
infinitesimal  (linear  problem),  but  also  not  excessively  large  (change  in  geometry,  very  large 
deflection  problem); 

•  Material  is  linear  elastic,  isotropic. 

1.2.2  Displacement  Functions 

The  finite  element  method  assumes  that  the  displacement  solution  is  a  node  displacement  vector 
fV .  For  a  plate  structure,  the  displacement  vector  consists  of  flexure  and  membrane  displacement 
vectors  Wj.  and  W„ ,  i.e. 


\wA 

(1-1) 

Relatively,  the  element  displacement  vector  can  be  expressed  as 

T 

II 

(1-2) 

The  displacements  at  any  point  on  the  element  could  be  expressed  in  terms 
displacements  as 

of  element  nodal 

w  =  C^Wj. 

(1-3-a) 

«  =  Cu^n. 

(1-3-b) 

(1-3-c) 

where,  the  ,  and  and  C„  are  the  row  vectors  of  interpolation  fimctions. 


1.2.3  Nonlinear  Strain-Displacement  Relation  (von  Karman) 
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Considering  small  inplane  strain  and  large  lateral  displacement,  the  total  strain  vector  is  given  by 
c  =  e  +  z  (1-4) 

where  the  membrane  strain  vector  e  consists  of  two  parts: 

e=e„+eg  (1-5) 

The  linear  membrane  strain  vector  is  related  to  the  displacements  as 


(1-6) 


Here  a  comma  represents  differentiation.  The  nonlinear  stretching  strain  vector,  ,  induced  from 
large  lateral  deflection  by  the  von  Karman  strain-displacement  relation  is  given  as 


The  vector  of  bending  curvatures  ?  in  equation  (1-4)  is  expressed  as 


(1-8) 


By  using  the  finite  element  displacement  functions,  equation  (1-3),  the  membrane  strain  and 
curvature  vectors  can  be  expressed  in  terms  of  the  element  nodal  displacements.  The  linear  membrane 
strains  from  equations  (1-3-b),  (1-3-c)  and  (1-6)  are  given  by 


The  nonlinear  membrane  strains  from  equations  (1-3-a)  and  (1-7)  are  expressed  as 
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0  ' 

r  -\ 

\Wx 

1  .X 

0 

yy 

w 

\x 

L  >yj 

=-?e=-? 
2  2 


—C 
dx  " 

— C 

ay 


Wf=-?C,Wf 


where  the  slope  matrix  and  vector  ^e  given  by 
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0  ' 

fX 

0 

Wy 

,  d  =  \  ’M 

>y 

w 

_,y 

w 

and  the  curvatures  from  equations  (1-3-a)  and  (1-8)  are  expressed  as 


-ilc 

dx 

?  =1 — 

ay'  " 

ilc 

w 

axy 


Wy  = 


1.2.4  Relations  of  Stress  and  Strain 

The  general  stress-strain  relation  for  a  plane  stress  (ct^  =  0 )  is  given  by 


s  = 


'V 


\=Ee 


where  the  strain  vector  is  given  by  equation  (1-4)  and  the  elastic  coefficient  matrix  is 

0 


E  = 


1-v' 


1  V 
V  1 
0  0 


0 

1-v 

2  J 


The  membrane  or  inplane  force  vector  N  and  moment  vector  M  are  defined  as 

/  2 

{N,M)=  j  {l,z)adz  =  {Ae,  Dk) 

where  A  and  D  are  symmetirc  matrix  of  material  properties;  h  is  the  thickness  of  the  plate, 
isotropic  plate  of  uniform  thickness,  they  are  given  by 


(1-10) 


(1-11) 


(1-12) 


(1-13) 


(1-14) 


(1-15) 
For  an 
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1  V  0 


1  V  0 


A  =  -^v  1  0  ,  ,  V  1  0 

0  0  ’  0  0 

2  . 


1.2.5  Equations  of  Motion 

According  to  the  principle  of  virtual  work,  the  total  work  done  by  internal  and  external  forces 
(including  inertia  force)  on  an  infinitesimal  virtual  displacement  is  zero,  that  is, 

SfF  =  dfVi„,-SfV^  =0  (1-17) 


where,  the  internal  and  external  work  are  given  by 

(1-18) 

Substituting  all  the  corresponding  equations  into  equations  (1-17)  and  (1-18),  the  element  equilibrium 
can  be  reached  as  (not  including  damping  here) 

Wff  0  n-iQi 


where  and  are  the  mass  matrices  of  flexure  and  membrane,  respectively.  They  are  defined 


Psh\plC„ds ,  m„„,  =  pAflCuvds 

*/.  =  ‘i.  =  Ufif’AC.ds .  t^=kl,=UciA?  C. 


kff,  k'ff  and  k^j-  are  the  linear,  the  first  order  nonlinear  and  the  second  order  nonlinear  parts  of  the 
stif&iess  matrix.  They  are  defined  as 

=  [cjDC^ds ,  4  =  ^[C^NSeds  ,  k\  =  ^ J Cj A?  C,ds  (1-22) 

and  mid-surface  force  matrix 


TV..  = 


N  N 

mx  nvcy 

N  N 

’  mxy  ^  my 


is  constructed  from  ’ 
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1.3  DKT  Triangular  Element 

In  the  following,  the  DKT  element  shown  in  Figure  1.1  is  adopted  from  Batoz  et  af  , 
Jeyachandrabose  et  af^  and  Xue“  without  modification.  DKT  element  defines  element  shape 
fimctions  due  to  slopes  as^® 


x,u 

Figure  1.1  DKT  triangular  element 


dx  oy 

where  ^  and  rj  are  area  coordinates  and  the  displacement  vector  is  expressed  as 

Wy  =  W,,  W^^,  Wj,  W^2>  '^3’  '^x3’  ^ 

The  nine  components  of  shape  fimction  vectors,  and  ,  are  given  by 

H^,=-\.5  a,N,-a,N, 


(1-26) 

(1-27) 

(1-28-a) 

(1-28-b) 
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0-28-c) 

=  0-28-d) 

H,2=-bsNs+kN,  a-28-e) 

Hy,  =  N,-e,N,-e,N,  0-28-f) 


The  functions  and  are  obtained  from  the  above  expressions  by 

replacing  A^,  by  N2  and  indices  6  and  5  by  4  and  6,  respectively.  The  functions  ,  H^g ,  , 

Hyg ,  Hyg  and  Hyg  are,  respectively,  obtained  by  replacing  by  Ng  and  indices  6  and  5  by  5  and 


4.  The  shape  functions  are  given  by 

=  2(1  - 77  ^  -t?  j  (1-29-a) 

N 2=  ^2^-1 

7^3=77277-1  (1-29-c) 

Ng=4^T]  (1-29-d) 

7\7j  =  477  1  -  4  -77  (1-29-e) 

TVj  =4<§  l-i§ -77  (1-29-f) 


and  the  coefficients  are  given  by 
~  ~  ^i)l  ^ij 


1  2  1  2 


V 


ih4+yl 


where  k=4,5,6  for  the  sides  i  =  23,  31, 12  respectively,  and 
Xy=x,-Xj,  yij  =  yi-yj 


(l-30-a) 

(1-30-b) 


(1-30-c) 

(1-30-d) 

(1-30-e) 

(1-30-f) 

(1-31) 
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According  to  equation  (1-10),  the  slope  transformation  matrix  Cg  can  be  found  as 


Q  = 


H. 


Cr  = 


1 

2A 


The  curvature  transformation  matrix  C ^  can  be  derived  as 

—  X^\H y  ^  “  X\2Hy  J^ 

—  X^yH^^  —  X^2^x,7]  y2\^y,^  y vfl y.-q 

where  A  is  the  area  of  an  element 

2  A  =  ^3i>’|2  “  ^12^31 

The  matrix  and  „  are  given  by 


<1  = 


Pel-2^  +ri  p^-Pi 

4  -  6(^  +  ■n)-r^il  -  2^ )+ 77(^5  +  rj 

96(1-2^)-T7{^5+96) 

-P6{i-2^)  +  V{P4  +  P6) 

2  -  6^  -  rg(l  -  2<^)-7](r4  -  rj 

II 

^6i^-24)-viq6-<i4) 

-niPs+P*) 

V{rs-r4) 

1 - 

1 

L4* 

- 

-Psl-2ri P,-Ps 

“ 

4  -  6((^  -t-  77 )  -  rj  (1  -  27] ) + 4  (r^  +  rj 

95(1-277)-^(^5+^J 

^{P4  +  P6) 

<^(^6-^4) 

.  ^4  = 

^(94 -?6) 

P5{^-2n)-^{P4+P5) 

2-6ri-r,il-2ri)-^{r^-r,) 

qs(l-2v)+^(q4-qs) 

where, 


t^\-2^  +r]  t^-t^ 

96(1-2^)-^(^5+?6) 

1  +  ^6(1-24)-17(''5  +  ''6) 
-U{\-2^)+r]{U  +  U) 
9g(l-2^)-|-7?(^4-96) 
-l-l-rg(l-2^)-l-T7(r4-rg) 

-niu+u) 

n(f4-f5) 

—  ^5 1  —  2tj  — 
qi{\-2ri)-^{q,+q^) 
l  +  r5(l-277)-(^(r5+rJ 

^(^4-96) 

^(^4-^6) 

^5(1  —  27])—  <^(^4  +  tj) 
?5(i-2n)+^(?4-95) 


Pk=-^xjll,  t,=-6yjll,  q,=^XyyJll,  r,  =  3yllll 

A:  =4,  5,  6  for  /  =23,  31,  12  respectively. 


(1-32) 


(1-33) 


(1-34) 


(1-35-a) 


(1-35-b) 


(1-36) 
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1.4  Statically  Nonlinear  Problems  of  Large  Deflection 

1.4.1  Scheme  for  solving  nonlinear  static  equations 

The  structural  static  equation  is  given  by 
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Equation  (1-40)  is  equivalent  to  tiie  following  two  equations 

KJtV.  -K^W, 


Hence,  the  iterative  equations  are  given  by 

Wp  =  Ff-K%W^;-'^  fVp'^ 
Wf  =  RfWf  +  (l  - 

.  WP  =  RJVP  +  {\-RjWp'^ 


0=1,2,-) 


and 


wp=0,  wp  =  0 


where,  Rj^  and  R„  are  relaxation  factors  for  the  flexure  and  membrane. 

The  scheme  for  solving  the  nonlinear  static  problem  is  listed  as  below. 

1.  Select  constants  for  iteration; 

2.  Formulate  the  load  vectors  Fj-  and  F^ ; 


3.  Initialize  the  displacement  vectors  fVj-  and  W„ ; 

4.  Calculate  the  stiffness  matrix  of  membrane  and  then  factorize  it; 

5.  For  ith  iteration: 

5.1.  Calculate  the  stiffhess  matrix  and  general  force  vector 

5.2  Solve  the  first  equation  of  equation  (1-42)  to  obtain  ; 

5.3  Check  convergence  for  flexure; 

5.4  Construct  new  flexure  displacement  vector  fV^  by  using  Rj^; 

5.5  Calculate  the  general  force  vector  ; 

5.6  Solve  the  third  equation  of  equation  (1-42)  to  obtain  WP  ; 

5.7  Check  convergence  for  membrane; 

5.8  Constmct  new  membrane  displacement  vector  WP  by  using  R„ ; 

5.9  If  the  displacements  of  flexine  or  membrane  are  not  convergent,  go  back  to  5.1; 


(1-40) 


(1-41-a) 

(141-b) 


(1-42) 


(1-43) 
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6.  Ou^ut  results. 


1.4.2  Results  and  Discussions 
1.4.2.1  Accuracy  of  the  Program 
Example  I 

A  clamped  squared  plate  subjected  to  a  uniform  distributed  load  is  considered.  The  length  of 
the  sides  is  7.62m.  The  thickness  is  0.0762m.  Modulus  elasticity  and  possion  ratio  are 
E  =  2.0685  xlO"N/m^  and  0.316  respectively.  This  example  has  become  a  “standard”  used  by  many 
authors^^'^'  to  check  the  GNL  (Geometric  Nonlinear)  thin  plate  formulations.  The  analytical  thin  plate 
solution  is  given  by  Levy^^  who  solved  von  Karman’s  equations  using  a  double  Fouries  series.  This 
solution  is  quoted  as  having  a  possible  error  of  less  than  2%^®.  The  plate  is  divided  into  200  triangular 
elements  and  121  nodes.  The  error  tolerances  are  10"^  and  Rj^  =  0.7 .  The  results  are  shown  in  Table 

1.1.  Q  =  (qa^)l(Et*)  is  a  non-dimensional  uniformed  distributed  load.  The  third  colunrn,  shown  in 

Table  1.1,  was  obtained  from  Wood^®  by  using  the  Irons-Razzaque  non-conforming  triangular  thin 
plate  element.  The  results  in  Columns  4  and  5  are  based  on  Mindlin’s  plate  theory.  Obviously,  the 
results  of  present  method  are  acceptable. 


Table  1.1 .  Comparison  of  central  deflections  from  different  methods  (clamped  plate) 


Load  (0) 

Exact 

Wood 

Pica,  et  al 

Rao,  et  al 

Present 

17.79 

0.237 

0.2387 

0.2351 

0.2346 

0.2383 

38.3 

0.471 

0.4717 

0.4673 

0.4660 

0.4636 

63.4 

0.695 

0.6916 

0.6887 

0.6862 

0.6791 

95 

0.912 

0.9008 

0.9003 

0.8987 

0.8721 

134.9 

1.121 

1.1025 

1.1041 

1.1025 

1.0677 

184 

1.323 

1.2961 

1.2990 

1.2979 

1.2595 

245 

1.521 

1.4879 

1.4913 

1.4890 

1.4539 

318 

1.714 

1.6744 

1.6774 

1.6750 

1.6475 

402 

1.902 

1.8529 

1.8682 

1.8526 

1.8365 

Example  11 

The  program  is  checked  again  by  using  a  simply  supported  square  plate  with  immovable  inplane 
edges.  The  other  parameters  are  identical  to  the  above  example.  The  results  are  listed  in  Table  1.2.  The 
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results  in  the  second  column  were  obtained  by  Rushton^^  who  used  a  finite  dilference  dynamic  relation 
approach  based  on  thin  Kirchhoff  plate.  There  is  a  little  difference  between  these  results. 


Table  1.2.  Comparison  of  central  deflections  fi-om  different  methods  (simply  supported) 


Load  fO) 

Rushton 

Pica,  et  al 

Present 

9.16 

0.335 

0.3478 

0.3296 

36.6 

0.818 

0.8184 

0.7625 

146.5 

1.47 

1.4655 

1.4168 

586 

2.40 

2.3927 

2.4520 

2344 

3.83 

3.8124 

4.1733 

9377 

6.07 

6.0521 

7.1272 

1.4.2.2  Effects  of  Rj-  on  the  convergence 
Example  I 

In  this  analysis,  a  simply  supported  square  plate  subjected  to  an  uniformed  distributed  load  is 
considered.  The  length  of  the  sides  is  10m.  The  thickness  is  0.1m.  Modulus  elasticity  and  possion  ratio 
are  E  =  2.0x10"  nW  and  0.3  respectively.  The  plate  is  divided  into  200  triangular  elements  and 
121  nodes.  The  iterations  and  central  deflections  for  different  are  listed  in  Table  1.3.  The  error 

tolerances  are  0.01.  Obviously,  R^  should  be  close  to  1  for  linear  problems,  while  it  is  much  smaller 
for  nonlinear  problems. 


Table  1.3.  Iterations  and  central  deflections  for  different  Rj^ 


Pressure 

(Pa) 

Deflection  of 

Linear 

/?^=1.0 

r=Q.i 

o 

ii 

Iter. 

Nonlinear 

Iter. 

Nonlinear 

Iter. 

Nonlinear 

200 

.440689e-03 

1 

.440673e-03 

5 

.440673e-03 

10 

.440673e-03 

2000 

.440689e-02 

1 

.439068e-03 

5 

.439089e-02 

9 

.439114e-02 

20000 

.440689e-01 

B 

.355715e-01 

2 

.356139e-01 

.356482e-01 

100000 

.220345e+00 

. 

15 

.896526e-01 

.900227e-01 

200000 

.440689e+00 

- 

_ 

- 

B 

.1225996+00 

2000000 

.440689e+01 

- 

- 

- 

- 

11 

.304909e+00 

Example  II 

The  example  II  in  section  1. 4.2.1  is  considered  again  to  check  the  convergence  of  the  iteration. 
The  iterations  for  different  cases  are  listed  in  Table  1.4.  The  results  show  that  it  is  very  difficult  to  find 
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a  single  optimal  Rj  for  different  load  cases.  For  higher  nonlinearity,  lower  7?^  is  needed  as  shown  in 
Table  1.4. 


Table  1.4.  Iterations  for  different  methods  (simply  supported  plate) 


Load  (Q) 

o 

II 

R=0.2 

/f/=0.05 

9.16 

12 

38 

169 

36.6 

11 

34 

148 

146.5 

16 

34 

150 

586 

31 

38 

159 

2344 

- 

52 

158 

9377 

- 

- 

201 

1.5  Dynamic  Nonlinear  Problems  of  Large  Deflection 


1.5.1  Scheme  for  solving  the  nonlinear  dynamic  equation 

The  classical  Newmark-  p  method^''  is  used  to  approximate  in  time.  For  convenience,  the  damping 
is  not  considered.  The  dynamic  equations  of  the  large  deflection  problem  without  damping  can  be 
written  in  matrix  form  as 


\Mff  0  -]\wA 

[  0 


K 


fm 


'^/L  P/ 

w„\  \f„ 


According  to  tiie  Newmark  method,  the  effective  stiffness  matrix  at  time  t  +  At  is 
Substituting  equation  (1-44)  into  equation  (1-45),  one  has 


/+A/ 

1 

> 

t+At 

7 

Km 

t+At 

+ 

0 

Kn,_ 

U 

L  0 

Because  Mj-j-,  and  are  constant  matrices,  we  have 


^JT~  "// 

t+At  1^  _ t-i-At  -w^ 

t+At  ^t  +  At  jcr 
t+At 

n  ^ 

mm 


(1^) 


(1-45) 


(1-46) 


(M7) 


20 


CHAPTER  1:  Nonlinear  Plate  Analysis  Using  Total  Lagrangian  Formulations 


According  to  the  Newraark  method,  the  effective  load  vector  is  given  by 
+  m{^‘W  +  +  a,'w) 

which  can  be  rewritten  in  partitioned  form  as 


l+Al  r  „ 

F. 


I+AI 


/u 


M 


0 


ff 

0  M. 


.Kia- 


•+a. 


Hence,  the  effective  load  vectors  of  flexure  and  membrane  are  given  by 

+  Mj^  a,  ‘Wj-  +  a^‘W^  +  a,‘W^ 

+  a,‘fV^+a,‘W„+a,‘W„ 

The  final  equations  to  be  solved  are  given  by 


t+Atr  ^ 


K 


ff 


^fm 


\_^mf  ^mm 


t+At 


w, 


t+it  r  « 


/u 


\fV„ 


According  to  the  assumption,  one  has  =  0 .  Therefore,  the  iterative  equations  for 

equations  are  given  by 


(1-48) 

(1-49) 

(l-50a) 

(l-50b) 

(1-51) 

effective 


t+At]^{i)  t+Atjyii-l)  t+Atjj^(i)^t+At p  ^t+Atp{0  t+Atjy(i-\)  t+AljyO^^) 


-h  (l-  i?y) 

^  mm  m  m  mf  \  f  f  J 


(0 

/ 


(/  =  U,-0 


(1-52) 


The  scheme  for  solving  the  nonlinear  dynamic  problem  is  listed  as  below. 

1.  Select  constants  for  Newmark-  p  method  and  the  static  iteration; 

2.  Calculate  integration  constants  for  Newmark-  p  method; 

3.  Initialize  structural  data  and  evaluate  index  vectors  for  diagonals  in  the  structural  matrices  of  flexure 
and  membrane; 

4.  Initialize  the  displacement,  velocity  and  acceleration  vectors; 


5.  Evaluate  and  factorize  the  stiffiiess  matrix  of  membrane  ; 


6.  Evaluate  the  mass  matrix  of  flexure  M  jy ; 

7.  Calculate  the  dynamic  responses  step  by  step  (Newmark-  p  loop): 
7.1  Formulate  external  force  vectors  Fy  and  ; 

12  Evaluate  effective  loads  at  time  ^  for  flexure; 
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7.3  Initialize  displacement  vectors; 

7.4  For  the  ith  iteration  of  nonlinear  effective  static  problem: 

7.4.1  Calculate  stif&iess  matrix  and  general  force  vector  of  flexure; 

7.4.2  Calculate  effective  stiffiiess  matrix  forflexiure; 

7.4.3  Solve  the  first  equation  of  equation  (1-52)  to  obtain  ; 

7.4.4  Check  convergence  for  flexure; 

7.4.5  Construct  new  flexure  displacement  vector  fV^'^  by  using  R^; 

7.4.6  Calculate  general  force  vector  j 

7.4.7  Solve  the  third  equation  of  equation  (1-52)  to  obtain  ; 

7.4.8  Check  convergence  for  membrane; 

7.4.9  Construct  new  membrane  displacement  vector  Wj'^  by  using  R„ ; 

7.4.10  If  the  displacements  of  flexure  or  membrane  are  not  convergent,  go  back  to  7.4.1; 

8.  Output  results. 

1.5.2  Numerical  Examples 
1.5.2.1  Example  I 

A  simply  supported  square  steel  plate  with  immovable  inplane  edges  is  considered.  The  length  of 
the  sides  is  1.0m.  The  thickness  is  0.01m.  Modulus  elasticity  and  possion  ratio  are  E=2.0E1  INW  and 
0.3  respectively.  The  plate  is  divided  into  200  triangular  elements  and  121  nodes.  Assume  a 
concentrated  load,  =10  Sin20  N,  acts  on  the  center  of  the  plate.  The  central  responses  for 
different  time  step  sizes  (0.0004,  0.0002,  and  0.0001)  and  error  tolerances  (10‘^,  lO"^,  and  10  *)  are 
drawn  in  Figure  1.2  respectively.  The  iterations  for  these  three  cases  are  5,  8,  and  12  respectively.  The 
results  show  that  the  error  tolerance  has  much  effect  on  the  accuracy  of  the  responses.  When  the  error 
tolerance  is  large  the  responses  change  with  the  time  step  size.  For  this  problem,  the  error  tolerance 

can  be  set  10"*. 

I.5.2.2.  Example  II 

The  plate  used  in  example  I  is  considered  again.  Assume  the  concentrated  load  is 
=  5  Sin  CO  kN.  The  central  responses  of  the  plate  for  ft)  =  3  rad/s,  ft)  =  100  rad/s  and  ft)  =  200rad/s 

are  drawn  in  Figure  1.3  respectively.  The  error  tolerance  is  2x10"^ .  A  and  B  denote  the  responses 
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obtained  from  die  nonlinear  and  linear  theories  respectively.  For  the  nonlinear  problems,  19  to  23,  13 
to  14,  and  13  to  14  iterations  are  required  in  the  three  cases  respectively. 

1.5.2  J  Example  III 

Another  simply  supported  square  aluminum  plate  with  immovable  inplane  edges^’  is  considered. 
The  length  of  the  sides  is  0.8m.  The  thickness  is  0.025m.  Modulus  elasticity  and  possion  ratio  are 
E=6.9elON/m^  and  0.33  respectively.  A  moving  load,  P=1000^Ar,  travels  at  constant  velocity 
v  =  153.7w/5  in  the  y-direction  along  the  centerline.  The  plate  is  divided  into  200  triangular 
elements  and  121  nodes.  Rj^  =  0.7  is  used  here.  The  responses  at  the  center  of  the  plate  for  linear  and 
nonlinear  are  given  in  Figure  1.4.  Again,  A  and  B  denote  the  responses  obtained  from  the  nonlinear 
and  linear  theories,  respectively.  In  this  figure,  the  error  tolerance  is  2x10”^  and  the  iterations  are  13 
to  14  for  the  nonlinear  case. 


1.6  Non-Dimensional  Form  of  FEM  Formulations 

The  finite  element  formulations  defined  above  are  dimensional.  As  we  know,  dimensionless 
parameters  are  usually  used  in  the  fluid  equations.  To  make  them  compatible,  these  FEM  formulations 
are  to  be  non-dimensionalized  with  respect  to  the  length  of  the  plate  and  other  parameters  from  fluid. 
The  dimensionless  dynamic  equations  of  equilibrium  using  finite  element  method  may  be  obtained  by 
non-dimensionlizing  either  the  partial  differential  equations  of  nonlinear  plate  or  the  integration  form 
of  the  finite  element  formulations.  The  details  of  the  former  will  be  shown  in  Section  5.4.  The  non- 
dimensional  form  of  the  latter  is  to  be  provided  in  the  following. 

Assume  the  dimensionless  parameters  are  given  by 


~u~v~w~x 
u=-,  v=-,  W  =  j,  X=j, 


t  =■ 


tv^ 


P 

PaVl 


(1-53) 


where  /  is  the  length  of  the  plate;  is  the  velocity  of  the  air  in  the  far  field;  denotes  the  mass 
density  of  air;  is  the  air  pressure  acting  on  the  plate.  With  this  assumption,  the  relation  of  the 
membrane  displacement  vectors  between  the  dimensional  form  and  dimensionless  form  is  given  by 

(>-«) 


23 


Computational  Mechanics  Laboratory,  University  of  Arkansas  at  Fayetteville 


Selvam  and  Qu 


For  the  triangular  element  used  above,  the  dimensionless  membrane  and  flexural  displacement  vectors 
are  defined  as 


=  “i  ^1  “2  %  “3  ^3  (1-55) 

Wj.  =  w^,  w^,  >v,2,  tVj,  w^3,  ^  (1-56) 

Clearly,  the  dimensionless  flexural  displacement  vector  does  not  have  the  simple  relation,  as  shown  in 
equation  (1-54),  with  its  dimensional  form.  To  simplify  the  derivation,  defined  the  flexural 
displacement  vector  as 


Wf=  w^,  ,  W2,  /w,2.  ^2.  M's ,  lw^3, 


(1-57) 


Now,  we  have  the  relation 

Wj-  =  Wfl  (1-58) 

Using  the  definition  in  equation  (1-57),  the  following  relations  of  the  interpolation  functions  or 
other  parameters  between  the  dimensional  form  and  the  dimensionless  form  may  be  obtained  from 
equations  (1-3),  (1-9),  (1-11),  (1-10),  and  (1-12)  respectively,  that  is. 


C„  =  C„,  C,=Q 
C.=Cjl 
?  -? 

Ce=Cjl 

c,  =  c,lP 


(1-59) 

(1-60) 

(1-61) 

(1-62) 

(1-63) 

(1-64) 


Introducing  equations  (1-59)  through  (1-64)  into  equations  (1-20),  (1-21),  (1-22),  and  (1-25)  leads  to 
"*//  =  ClCJs  =  p^hPin ,  /«„„  =  p^hP |  ClC^^ds  =  p,hPm„„  (1-65) 


^ _ L  f  DC  ds  = 


(1-67) 
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l‘^  =  D.[clDCJs  =  D.K 


(1-70) 


where. 


D  = 


1  V  0 
V  1  0 
0  0 


= 


Eh 


1-v 


2  ’ 


2  J 

D,= 


Eh^ 


^  12(1 -v') 


N  = 


N_  = 


N  N 

nvc  mxy 

N  N 

^  ’  mxy  my 


N„ 


N. 


my 


N. 


mxy 


=  DC„w„ 

'  m  m 


(1-71) 

(1-72) 

(1-73) 

(1-74) 


Using  equation  (1-53),  the  dimensionless  acceleration  and  force  vectors  are  given  by 

~  ..  Fi  ~ 

M7„=-y-W„ 

f  =  pS  =  p/li^rS=PaVll^f 


(1-75) 


(1-76) 


Remember,  the  vector  Wy  in  equation  (1-75)  is  corresponding  to  the  displacement  vector  in  equation 

(1-57).  Inserting  equations  (1-65)  through  (1-70),  (1-75),  and  (1-76)  into  equation  (1-19)  results  in 

/ 


Pshiv: 


m 


ff 

0  in 


D 


I 


fr®  0 
Kjy  V 

0  0 


+  DJ 


k'  +k^  k 

Kjjf  -t  Kjj-  «y„ 


JJKJ  [L] 


(1-77) 


Dividing  both  hand  sides  of  equation  (1-77)  by  p^hlVl  yields 

1  /  r 


m 


ff 

0  in 


0 


w, 


7  1 


^ff  ^ 
0  0 


^ffj  ^ff  ^fm 

k  ^ 


-■7 


’4=pJ^  (1-78) 

mj  \jm J 


in  which  the  dimensionless  dynamic  pressiure  A ,  mass  ratio  ,  and  internal  force  coefficient  t„,  are 
defined  as 


X-Biy  n-&L  ,  =11 

D,  -  '-‘-I 


y>‘, 


(1-79) 
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Equation  (1-78)  is  the  dynamic  equilibrium  of  an  element  in  non-dimensional  form. 

Because  the  flexural  displacement  vector  used  here  is  different  from  that  used  in  the  triangular 
element  in  Section  1.3,  the  functions  pertaining  to  this  flexural  displacement  should  be  different  from 
those  defined  in  Section  1.3.  The  former  may  be  obtained  simply  by  dividing  the  items  pertaining  to 
the  rotational  degrees  of  freedom  in  the  latter  by  the  length  /,  namely, 

•  Flexural  displacement  functions: 

Zi^(l  +  24  +  2i3)+2Ii44 

•^K'^31'^3  “^12) 

■^|(T31-^3  "TnA)"!"  ">’12) 

4(1  +  24 +21,) +  21,413 

4(^124  ~'^234)'^t4'^24('^12  “^23) 


4(Ti24  “3'234)‘*‘ 2  4A45(Ti2  T23) 

4(i+24  +  24)+2444 

4(^23^  —^14)  + ^444(^23  —•*31) 

4(^234  ~T3i4)'*"^ 444 (3*23  “T31) 


Slope  transformation  matrix: 


(1-80) 


(1-81) 


in  which  and  are  the  shape  fiinction  vectors.  Their  components  have  the  same  forms  as 
those  shown  in  equation  (1-28)  while  the  coefficients  4 .  and  are  given  by 


4 


Ck  = 


1  2  1  2 


(  1  2  1  2 

e*-  j^yij 


ni 


(l-82a) 

(l-82b) 

(1-82c) 


Curvature  transformation  matrix  has  the  same  form  of  that  in  equation  (1-33)  except  the 
coefficients  and  in  equation  (1-36)  which  are  given  by 


26 


CHAPTER  1:  Nonlinear  Plate  Analysis  Using  Total  Lagrangian  Formulations 


Finally,  the  dimensionless  forms  of  these  functions  may  be  obtained  by  using  the  relations  in 
equations  (1-59)  through  (1-64).  They  are  listed  in  the  following  for  reference. 

•  Flexural  displacement  functions: 

4(1  +  24  +  2Z3)  +  2Z,4Z3  1"^ 

A  (^3 1 A  ”  ^12  A  AAA  (^1  ”  ■^12 ) 

A  (^31 A  J^nA  )'^  2  AAA  0^31  “^^12) 

40  +  213 +2Zi)  +  2Ii4l3 

^  ^  ^  A(*^12A  “ ’^23A) +tAAA(*^12  “ -^23)  ►  (a) 

~723^)'^'^AA'^(Ti2  ~  yii) 

I^{1+2L,+2L2)  +  2I^I^L, 

^{^23^2  ~^iA)'^‘^A-^2-^(-*'23  "-*31) 

^  (^23^2  “  >'3 1 A  )  +  ~  (>"23  —  3^3 1 ) 

•  Membrane  displacement  functions: 

~  fZ,,  0  4  0  Lj  Ol 

C  =  (b) 

[O  Z,  0  4  0  4 

•  Membrane  strain  coefficients: 

[>'23  ?3i  ?i2  0  0  01 


>'23 

?3I 

?12 

0 

0 

1 - 

0 

0 

0 

0 

X32 

Xj3 

(C) 

X32 

X,3 

^21 

J23 

731 

7i2. 

Slope  transformation  matrix: 

~  [H.l 


where  and  are  the  dimensionless  shape  function  vectors.  Their  components  have  the 

same  forms  as  those  shown  in  equation  (1-28)  while  the  all  the  variables  on  the  right  hand  sides  of 
equations  (1-30)  should  be  their  non-dimensional  forms. 

•  Curvature  transformation: 
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The  dimensionless  derivatives  of  the  shape  function  vectors  have  the  same  forms  as  those 
expressed  in  equation  (1-35)  except  the  definitions  in  equation  (1-36).  In  the  new  definitions,  the 
dimensional  variables  in  equations  (1-36)  are  replaced  by  the  dimensionless  variables  shown  in 
equation  (1-53). 


1.7  Summary 

The  large  deflection  of  thin  plates,  which  considers  the  coupling  effect  between  the  flexure  and  the 
in-plane  deformation,  was  investigated  by  using  finite  element  method.  Based  on  the  von-Karman 
nonlinear  strain-displacement  relation,  the  dynamic  equations  of  large  deflection  plates  were  derived  at 
first.  The  Formulation  of  Discrete  Kirchhoff  Theory  (DKT)  for  triangular  elements  was  then  presented 
for  the  discretization  of  the  thin  plate.  An  iterative  method  between  the  equations  of  the  flexure  and 
the  membrane  was  used  to  solve  the  static  nonlinear  problem.  To  make  the  method  efficient,  a 
relaxation  factor  was  used  in  this  method.  The  Newmark-  p  method  was  applied  to  solve  the  dynamic 

equations  in  time  domain.  Several  numerical  examples  were  included  to  demonstrate  the  efficiency  of 
the  method  for  the  nonlinear  dynamic  problem.  The  following  conclusions  can  be  drawn  from  the 
theory  and  the  results  of  the  examples: 

1.  The  results  from  the  present  program  are  acceptable  by  comparing  with  other  typical  results. 

2.  The  relaxation  factor  is  very  important  for  the  convergence  of  the  iterative  method.  If  improper 
relaxation  factor  is  selected,  the  iteration  may  not  converge. 

3.  Rj  should  be  close  to  1  for  linear  problems,  while  it  is  close  to  zero  for  nonlinear  problems. 

4.  The  error  tolerance  has  much  effect  on  the  accuracy  of  the  responses.  When  the  error  tolerance  is 
large  the  responses  change  with  the  time  step  size. 
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CHAPTER  2 


NONLEVEAR  PLATE  ANALYSIS  USING 
COROT  ATIONAL  FORMULATIONS 


2.1  Introduction 

2.2  Nonlinear  Static  Analysis  of  Plates 

2.3  Nonlinear  Dynamic  Analysis  of  Plates 

2.4  Summary 

2.5  Appendix 


2.1  Introduction 


Geometric  non-linearity  is  increasingly  being  considered  in  structural  analysis  especially  for  those 
structures  that  undergo  large  displacement  with  very  little  actual  staining  occurring,  such  as  beams  and 
plates.  There  are  many  works  on  the  way  in  which  the  geometrically  nonlinear  problems  may  be 
formulated.  Three  main  methods  are  currently  employed,  the  Total  Lagrangian  (TL)  formulation,  the 
Updated  Lagrangian  (UL)  formulation,  and  the  Co-Rotational  (CR)  formulation.  Each  of  these 
formulation  has  its  inherent  drawbacks,  merits  and  advocates.  The  TL  formulation  in  which  the 
original  undeformed  configuration  is  taken  as  the  reference  configuration  is  only  really  effective  for 
relatively  small  rotation.  The  UL  formulation  where  the  configuration  at  the  end  of  the  last  iteration  is 
used  as  the  reference  configuration  can  be  path  dependent  and  can  also  show  very  slow  convergence. 
The  CR  formulation  which  makes  use  of  both  the  original  configuration  and  the  current  deformed 
configuration  to  formulate  the  system  matrices,  although  restricted  to  small  rotation  between 
iterations*’^,  exhibits  good  rate  of  convergence  and  thereby  reducing  the  computation  time.  In  this 
work,  the  CR  formulation  for  static  and  dynamic  analysis  of  plate  structures  using  DKT  triangular 
elements  is  presented.  A  number  of  numerical  examples  are  presented  that  fully  test  the  capabilities  of 
the  formulation.  The  CR  formulation  code  developed  in  this  study  can  be  used  for  more  understanding 
and  realistic  assessment  of  the  panel  flutter  problems. 
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2.2  Nonlinear  Static  Analysis  of  Plates 

2.2.1  Elements  in  the  Analysis 

Large  deflection  plate  problem  generally  involves  the  coupling  of  in-plane  and  out-of-plane 
effects.  In-plane  deformation  associated  with  the  constant  strain  triangle  with  six  degrees  of  freedom, 
two  at  each  node,  and  out-of-plane  deformation  associated  with  the  DKT  trinagle  witii  nine  degrees  of 
freedom,  three  at  each  node,  are  shown  in  Figure  2.1.  Local  displacement  has,  in  total,  fifteen  degrees 
of  freedom,  five  at  each  node,  of  the  triangular  element.  Local  and  global  degrees  of  freedom  are 
represented  by  D  and  d,  respectively.  Simbolically,  they  are  given  by 


Figme  2.1  Degrees  of  freedom  (in-plane  and  out-of-plane) 
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D  =  Ui  V,  Wi  y/,,  Wyx  U2  V2  W2  y/,,  y/^^  U3  V3  W3  y/,,  y/^, '' 

and 

d  =  U,  V,  W,  0,,  0y,  U3  V2  W2  9^  0y2  U3  V3  W3  9^  0y3  ^ 
T  T 

For  plane  stress:  [u]  =  u,  Uj  U3  ,[v]=  v,  V2  V3  , 

T 

For  bending:  [w]=  w,  0^,  0^,  Wj  0*2  ^>2  Wj  0^3  0^3  , 


Where  0^i 


and  0y,  =-tan 


and  so  on. 


(2.1) 

(2.2) 


2.2.2  General  Procedures  of  Large  Deflection  Analysis  of  Plates 

1.  Generate  the  element  stiffness  matrix  Kp  (6x6)  for  plane  stress  condition  (appendix-2)  and 

K  (9x9)  for  plate  bending  (Appendix-3)  and  assemble  them  to  find  the  element  stiffness 
matrix  (15x15)  in  local  co-ordinates.  It  is  obvious  that  the  element  stiffness  matrix  in 

global  coordinates  would  be  the  same  as  ; 

2.  Assemble  all  the  element  stiffiiess  matrices  to  obtain  global  structure  stiffness  matrix  K  .  The 
equilibrium  equation  in  global  coordinates  is  [K  d  =  r  where  d  is  the  vector  of  nodal 
displacement  components  and  [r  is  the  vector  of  extrenally  applied  nodal  forces.  Solve  the 
equilibrium  equation  to  obtain  the  global  displacement  vector  d  ; 

3.  Compute  the  element  stiffness  matrix  in  its  current  local  configuration  and  establish  the 

transformation  matrix  [T  (15x15)  as  derived  in  Appendix-1; 

4.  Determine  the  local  displacement  vector  D  from  global  displacement  vector  d  .  Local 

displacements  are  defined  as  the  displacement  of  the  nodal  points  of  the  deformed  element  from 
the  nodes  of  an  undeformed  reference  element  (l  -2  ~3  )  located  in  the  plane  as 

shown  in  Figure  2.2; 

5.  Compute  the  local  resisting  force  Q  =  [D]; 
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6.  Transform  the  resisting  force  Q  and  stiffness  matrix  to  global  orientation  by  [T]^  Q  and 

[T]^  K.  [T]  ,  respectively,  and  assemble  them  to  obtain  q  and  K  in  global  coordinates; 

7.  Compute  the  increments  in  global  displacement  Ad  by  solving  the  incremental  equilibrium 
equation  [K  Ad  =  [ Ar  ,  where  [ Ar  =  [r  -  q  ; 

8.  Determine  the  new  estimate  of  global  displacement  [d  =[d  +[Ad  and  repeat  steps  3 
through  7  until  Ar  becomes  arbitrarily  small; 

9.  Apply  a  new  load  increment  and  follow  the  steps  3  through  8  until  equilibrium  is  established  for 
the  new  load  increment. 

2.2.3  Coordinate  System 

Coordinate  systems  before  and  after  deformation  of  the  triangular  element  used  in  the  analysis  are 
shown  in  Figure  2.2. 

1.  x,y,z  :  the  global  coordinate  system  before  deformation.  Plate  element  initially  lies  completely 

on  x-y  plane.  Hence,  Xj  ,yi  is  the  global  coordinates  of  any  point  on  the  element  before 

deformation.  The  global  coordinates  of  nodes  1,  2  and  3  are  X|,y]  ,  X2,y2  and  X3,y3  , 
respectively,  with  z-coordinate  of  each  node  being  zero. 

2.  X,  Y  :  Local  coordinates  of  any  point  on  the  element  before  deformation.  The  local  coordinates 

of  nodes  1, 2  and  3  are  X,,Y,  ,  X2,Y2  and  X3,Y3  ,  respectively. 

Xi,Yi  =  x.,y,  -  x,,y, 

3.  x‘*,y‘^,z‘*  :  Global  coordinates  of  any  point  on  the  element  after  deformation.  The  global 
coordinates  of  nodes  1,  2  and  i  of  the  deformed  element  are  (xj,y^,z^),  {^2’y 

X3,y3,Z3  , respectively. 

xf,yf,zf  =(xi,yi,o)+  Ui,Vi,Wi 

4.  :  Local  coordinates  of  any  point  on  the  element  after  deformation,  The  local 

coordinates  of  nodes  1, 2,  and  3  are  ,  X2,  Y2  ,Z2  and  X3,  Y^  ,Z3  ,  respectively. 
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5.  X“ ,  Y“  :  Local  coordinates  of  any  point  on  the  undeformed  reference  element  whose  first  node 

coincides  with  the  first  node  of  the  deformed  element.  Z“  -coordinate  of  the  undeforemed 
reference  element  (r'-2"-3'')  is  zero  as  it  completely  lies  on  X‘-Y‘‘  plane. 


X 


Figure  2.2  Coordinate  systems  before  and  after  deformation 


2.2.4  Determination  of  local  displacement  D  from  global  displacement  d 

1.  Compute  global  coordinates  after  deformation  (xf,yf,zf),  and  (x^.y^.z^)  from 

the  global  displacement  vector 

d  =  u,  V|  w,  0,1  0y,  Uj  V2  W2  6;^  0y2  Uj  Vj  Wj  6^  6y3  ^  obtained  in  step-2  of 
Art.  1.2  and  the  input  global  coordinates  x,y,z  before  deformation.  The  relationship  stands  for 

xf.yf.zf  =(xi,yi,0)+  Uj.Vi.Wj 

2.  Rearrange  the  displacement  vector  d  as 
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where,  tan  J  =  -tan~*(0xi) 

X  dy 


X"  = 


ax“ 

ax'* 


1  + 


ay** 


1+ 


ax” 


V 


and 


1  +  - 


aVj 

ay" 

aVi 


aUj  aVj 

ay“  ax“ 


y“  = 


ay“ 

ax'* 


ay” 


1+ 


ax“ 


V 


1+ 


aVj 

ay“ 


A 


aUj  aVj 

ay"  ax" 


Similar  values  for  X y  =  andyy  =  can  be  computed  as 


X"  =■ 


ax" 

ay'’ 


aUj 

ay" 


ax" 


V 


1+ 


A 


1+ 


ay" 

aui 


aUj  aVj 

ay"  ax" 


y^  = 


ay" 

ay^ 


ax" 


1+ 


aUj 

ax'* 


V 


3Y“ 


A 


aUj  aVj 

ay**  ax'* 


Now,Uxe  = 


auj 


ax" 


ax"  ^ 

au: 

<ay»  > 

w  3  Vi 

fax""! 

aVj 

<ay"> 

ax" 

V  J 

+ - !- 

ay" 

[ax-J 

’  ax" 

[ax'J 

+ 

ay" 

[ax'J 

Uye  = 


aUj 


ay" 


^ay" ' 

au: 

fax"  ] 

w  9Vi 

"ay"^ 

aVj 

rax""! 

ay" 

j 

+ 

ax" 

ay" 

^  j 

’  9Y" 

ay" 

v  J 

"^ax" 

[3Y-J 

can  easily  be  computed. 


^  Uj  auj  aVj 

From  equations  A2.10  and  A2,l  1  (Appendix-2),  it  can  be  seen  that  the  values  of  , 

dX  I  X 


y" 


are  all  constants  and  they  are. 


aUj 


x^  =  “2  =^[(y2-y3)Ui+(y3-yi)U2+(yi  -y2)U3] 
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■^  =  03  =  ^  [(xj  -  y 2  )U,  +  (x  1  -  X3  )U2  +  {x2  -  X  JU 3 ] 

=  P2  =  ^ [(y2  -  Ys  )Vi  +  (ys  -  Yi  )V2  +  (Yi  -  Y2  )V3  ] 

=  p3  =  ^  [(^3  -  Y2  )Vl  +  (xi  -X3  )V2  +  (X2  -Xi  )V3  ] 

T  T 

6.  Rearrange  the  to  comer  order  to  have  D  . 

2.2.5  Results  and  Discussion 

Example  1 

A  clamped  square  plate  subjected  to  a  imiformly  distributed  load  q  is  considered.  The  plate  is  of 
7620  mm  x  7620  mm  having  a  thickness  of  76.2  mm.  Modulus  of  elasticity  (E)  and  Poisson’s  ratio  (v) 
are  2.0685  xlO’N/mm^  and  0.316,  respectively.  This  example  has  also  been  used  by  Selvam  and  Qu’. 
The  analytical  thin  plate  solution  is  given  by  Levy^  who  solved  von  Karman’s  equations  using  double 
fourier  series. 

The  plate  is  divided  into  288  triangular  elements.  The  niunber  of  nodes  was  169.  The  error 
tolerance  was  10■^  The  results  are  shown  in  Table  2.1  where  Q  =  qa'‘/Et'*  is  a  non-dimensional 
uniformely  distributed  load.  It  can  be  observed  that  the  present  analysis  estimates  the  results  well. 

Example  2 

The  program  is  checked  by  using  a  simply  supported  square  plate  with  immovable  inplane  edges. 
The  other  parameters  are  identical  to  Example  1.  The  results  are  listed  in  Table  2.2.  The  program  was 
checked  for  linear  analysis  using  a  concentrated  load  P  equal  to  lO’  N  at  the  center  of  the  plate.  The 

result  is  checked  against  the  standard  solution  (5^^^  =  a  Pa^/D  =7.9495  mm)  as  derived  by 
Timoshenko^.  The  maximum  deflection  was  found  to  be  7.9324  mm  with  an  estimated  error  of  0.2%. 
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Table  2.1  Comparison  of  the  ratios  of  central  deflections  to  thickness 


for  different  methods  (clamped  plate) 


Q 

Exact 

Wood 

Pica,  et  al 

D™.x/t 

Rao,  at  al 

Selvam,  et  al 

Present 

%  Error 

17.79 

0.2383 

^^9 

0.4673 

0.4660 

0.4636 

0.4748 

63.40 

0.695 

0.6916 

0.6887 

0.6862 

0.6791 

0.6990 

0.576 

95.00 

0.912 

0.9008 

0.9003 

0.8987 

0.8721 

0.9131 

0.121 

134.90 

1.121 

1.1025 

1.1041 

1.1025 

1.0677 

1.1193 

0.152 

184.00 

1.323 

1.2961 

1.2990 

1.2979 

1.2595 

1.3168 

0.469 

245.00 

1.521 

1.4879 

1.4913 

1.4890 

1.4539 

1.5106 

0.684 

318.00 

1.714 

1.6744 

1.6774 

1.6750 

1.6475 

1.7006 

0.782 

402.00 

1.902 

1.8529 

1.8682 

1.8526 

1.8365 

1.8802 

1.146 

Table  2.2  Comparison  of  the  ratios  of  central  deflections  to  thickness 


for  different  methods  (Simply  supported) 


Q 

Russton 

D 

Pica,  et  al 

maxA 

Selvam  et  al 

Present 

9.16 

wiu«ii| 

0.3296 

0.3424 

36.60 

0.7625 

0.8073 

1.4655 

1.4168 

1.4471 

2.3927 

2344.00 

3.8300 

3.8124 

4.1733 

3.7836 

9377.00 

6.0700 

6.0521 

7.1272 

6.0273 

The  computation  in  each  loading  was  obtained  by  using  a  single  load  increment  i.e.  the  load  for 
which  the  solution  is  sought  is  applied  in  full  at  a  time  and  the  iteration  was  performed  within  this  load 
until  the  difference  of  the  applied  full  load  and  the  resisting  forces  becomes  niinimum(of  the  order  of 
10"^).  The  computation  can  also  be  performed  by  using  as  many  load  steps  as  desired  without  using 
the  full  load  at  a  time.  But  the  number  of  iterations  is  much  less  for  the  former  case.  As  we  are 
interested  in  the  central  deflection  of  the  plate,  the  computaion  was  performed  with  the  error  tolerence 
based  on  the  difference  of  the  central  deflection  in  consecutive  two  iterations.  It  was  observed  that  the 
solution  begins  to  diverge  at  a  load  equal  to  and  greater  than  Q=586  (Dn,a)/t=2.314)  for  which  the 
problem  becomes  highly  non-linear.  For  convergence  of  the  solution  for  higher  loads,  a  relaxation 
factor  (Rf)  of  0.3  was  used.  After  a  few  trial  computations  it  was  found  that  the  result  begins  to  diverge 


45 


Computational  Mechanics  Laboratory,  University  of  Arkansas  at  Fayetteville  Selvam,  Roy,  and  Qu 

for  a  load  when  Dmai/t  becomes  approximately  equal  to  1.90  and  an  introduction  of  relaxation  factor  is 
then  needed  for  the  convergence  of  the  solution. 

Example  3 

In  this  case,  a  simply  supported  square  plate  of  size  and  thickness  equal  to  16  inch  and  0.1  inch, 
respectively  was  used.  In  the  computation,  the  modulus  of  elasticity  (E)  and  Poisson’s  ratio  (v)  were 
taken  as  30x10*  psi,  0.316,  respectively.  This  example  was  used  by  Murray^.  The  central  deflection 
was  carried  out  for  q=15  psi  in  5  increments  and  the  results  are  shown  in  Table  2.3  along  with  the 
results  read  from  the  graph  produced  by  Murray.  For  higher  loads  the  deflection  was  computed  and  the 
results  are  shown  in  Table  2.4. 


Table  2.3  Comparison  of  the  ratios  of  central  deflections 
to  thickness  (Simply  supported) 


q,psi 

Analysis 

D^ax/t 

Murray  (read  from  graph) 

3 

1.0499 

1.05 

6 

1.3868 

1.38 

9 

1.6122 

1.60 

12 

1.7875 

1.78 

15 

1.9335 

1.92 

Table  2.4  Central  deflection  of  the  plate 


q,psi 

20 

40 

60 

80 

100 

120 

140 

Dn,a=/t 

2.1362 

2.7044 

3.0990 

3.4129 

3.6774 

3.9088 

4.1152 

2.3  Nonlinear  Dynamic  Anaiysis  of  Plates 

2.3.1  Formulation  of  Nonlinear  Dynamic  Equations 

Newmark  implicit  time  integration  scheme  is  employed  to  calculate  the  nonlinear  dynamic 
response.  For  convenience,  the  damping  is  not  considered  in  the  present  calculation  although  the 
inclusion  of  damping  will  not  affect  the  formulation.  The  dynamic  equations  for  the  large  deflection 
problems  without  damping  can  be  written  in  matrix  form  as: 
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0 

'Db 

+ 

0  ■ 

d; 

- 1 

_ 1 

0 

Mp_ 

Pp. 

0 

K, 

Pp. 

(2.3) 


where  the  subscripts  B  and  P  represents  the  characteristics  associated  with  the  bending  and  in-plane 
nodal  displacements,  respectively.  Equation  (2.2)  can  be  written  in  general  form  as 

m  +  KD=F  (2.4) 

If  the  thickness  of  an  element  is  t  and  that  is  assumed  to  be  constant  widiin  the  element,  the  consistent 
mass  matrix  can  be  computed  from  the  following  equation; 

Mg  =  ptJjN^N  dxdy  (2.5) 

The  following  displacement  functions  are  used  to  develop  the  mass  matrix  for  flexiue’ 


Zr|  +  ^2  ”  A  ^2  ~  ^1^3 


1 


A  ^3  -  ^if'2^3 


-b. 


{^544 


V  “  '  '  ^ 

l“^3f  a  a 


2^3  I  ‘-3 

2  j-  I  r2 

1  ^  ^ 


^2  +  ^2^3  "^^2  A  ~  A  A  ~  A  A 


A  A  XA1A2A 

V  ^ 

^  .  1  ^ 

A^+~AAA 

\  ^  j 


\ 


\  ^2 

^  L 
1 


-c( 


A2A  ■^^1X2  A 
4l3+il,4l3 


Z,  +Z'I,  +L\L,  -X,Z?  -L,L. 


3^2 


b.  4%+-z.4Z3 


+i 


X1X2X3 


/  1  ^ 

A3X2  +  2  A1A2A3  J-C2 


(2.6) 

V  -  yj 

The  mass  matrix  for  flexure  has  been  computed  numerically  following  the  procedures  given  by 
Anderson®  and  that  for  membrane  can  be  determined  explicitly  and  is  given  by  Zienkiewic^.  The 
formulation  of  stiffness  matrix  has  been  described  in  the  static  part. 

Employing  Newmark’s  assumptions 

=  'Z)  +  [l-5  ‘D  +  S  (2.7) 

■^1 


‘*^D  =  ‘D  +  ‘DAt  + 


-a  \‘D+a‘*^‘D 


LV 


(2.8) 
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The  equilibrium  equation  at  time  t  +  At 

M  +  K  '*^D  =  (2-9) 

Expression  for  ‘*^‘D  can  be  obtained  by  solving  equation  (2.8)  for  '^^‘D  in  terms  of  the  unknown 
displacement  and  then  substituting  for  ‘*^‘D  into  equation  (2.7),  expression  for  the  velocity 


‘*^'D  can  be  obtained.  These  are  as  follow: 

D=^ao  ‘D  -a^'D- Oj ‘D  =  “  ‘Q 

D  =  <D  +  a^'D  +  ay  ‘*^‘D  +  ‘R  (2-11) 

where  '0  =  [ao 'D  +  Oj 'D  +  flj 'jD]  and  7?  =  [ 'D  +  a^ 'd].  Expression  for  the  unknown 
displacement  ‘*^‘D  can  be  obtained  by  substituting  equation  (2.10)  into  equation  (2.9)  as 

a^M  +  K  (2-12) 

Where  the  effective  load  vector,  ‘^‘^F  at  time  t  +  At  is  given  by 

'*^F  +  M  ao‘D  +  a2‘D  +  a^‘D=‘*^R+M  ‘Q  (2.13) 


2.3.2  iteration  schemes 

Using  implicit  time  integration,  the  equilibrium  of  the  system  at  time  t  +  At  is  considered.  This 
requires  in  nonlinear  analysis  that  iteration  be  performed.  Neglecting  the  effects  of  damping  matrix 
and  lining  the  modified  Newton-Raphson  iteration,  the  governing  equilibrium  equations  are: 

«+A(^n  _t+Mp _l*^l pn-\  (2.14) 


and 


/+A<  ^  (2.15) 

where  is  the  static  resisting  force.  Using  the  expression  of  acceleration  at  time  t  +  At, 

in  equation  (2.10),  equation  (2.15)  can  be  written  as 


and  the  velocity 

D" +‘R  (2-17) 

Substituting  equation  (2.16)  in  equation  (2.14),  the  final  expression  for  the  incremental  displacement 
AD”  can  be  obtained  as 
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2.3.3  Steps  of  Computing  Nonlinear  Dynamic  Response 

1.  Compute:  Stiffness  matrix,  K,  Mass  matrix,  M ,  Damping  matrix,  C; 

2.  Set  initial  conditions  at  time.  Displacement,  and  Velocity,  D^; 

3.  Compute  initial  accelerations,  Dq,  Dq  =M  ’[fO  —  CDq— KDq]; 

4.  Select:  Time  step.  At  and  Parameters,  a  and  y.  Calculate  Integration  constants  ao.  ai,  %,  34,  as, 

a6,  and  a?; 

5.  Form  effective  stiffness  matrix  K,  K  =  K-^  a„M  +  a,C  and  tutorize  it  K  =  LUL  ; 

6.  Calculate  effective  loads  at  time  t + At ,  F  =‘*^F  +  M  'Q  +  C‘S ,  where, 

'g  =  [a/D  +  02  ‘b  +  aj  ‘b]  and  ‘S  =  [a,  ‘D  ■\-a^‘b  +  ‘b] ; 

7.  Solve  for  displacement  at  time  t  +  At ,  from  {lUII  )  F ; 

8.  Calculate  acceleration  and  velocity  at  time  t  +  At ,  ‘*^^0  =  0^  ‘*^‘D—‘Q  and 
•*‘^b  =  a2‘*^b+‘R; 

9.  Now  for  the  n-th  iteration  at  this  instant  of  time,  t  +  At: 

(a)  Transform  the  displacement  '*^Dinto  local  co-ordinates  C^d),  establish  local  element 

stiffness  C^k)  and  element  mass  matrices.  Compute  local  balancing  elastic  forces, 

and  transform  and  assemble  it  to  global  axis  ).  Compute  the  global  mass 

matrix  from  and  global  stiffness  matrix  from 

(b)  Calculate  the  unbalanced  forces  at  time,  t  +  At  at  the  n-th  iteration  as 

f+A/  /+Af  t+Atj^n-\ 

(c)  Compute  the  effective  stiffriess  as  +aQ'*^'M" 

(d)  Solve  for  incremental  displacement AD"  AF  ” ; 

(e)  Update  displacement,  acceleration,  and  velocity  by 

+AD",‘*'^b''  =ao‘*^'D"-'-‘Q  +  aoAD'‘  and  ‘*^b"  =ao‘^^b"+‘R, 

respectively; 
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(f)  Continue  from  9(a)  to  9(e)  until  the  both  or  either  of  the  following  two  criteria  has  been  met 
(i)  an  error  limit  based  on  unbalanced  forces  i.e.  '  AF  /  F  <  specified  limit,  (ii)  an 
error  limit  based  on  displacement,  i.e.  '"’'“AD"/  <  specified  limit. 

2.3.4  Results 

Example  1: 

A  simply  supported  square  steel  plate  is  considered.  The  length  of  the  side  and  thickness  are 
respectively,  1.0  m  and  0.01  m.  The  modulus  of  elasticity  and  the  Poisson’s  ratio  are  2.0<10"N/m^ 
and  0.3,  respectively.  A  concentrated  load  P=50000sin  cot  N  is  assumed  to  be  acting  on  the  center  of 
the  plate  and  the  central  responses  were  computed  for  CD=  3  rad/s,  100  rad/s  and  200  rad/s. 

Figure  2.3  shows  the  central  responses  at  different  time  steps  with  different  error  levels  at  a 
forcing  fi^quency  of  100  rad/s.  Figure  2.4  shows  the  linear  and  nonlinear  central  responses  at  different 
frequencies.  At  low  frequencies  (3  rad/s)  the  nonlinear  responses  are  similar  to  results  presented  in  the 
report  of  Selvam  and  Qu,  but  at  higher  frequencies  (100  rad/s  and  200  rad/s)  the  response  does  not 
match  exactly  with  the  report  especially  at  later  times. 


(a) 
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(b) 


Figure  2.3  Central  responses  of  the  plate  at  03=3  rad/s  at  different  time  steps  with  an  error  level  of 

(a)  lO’  and  (b)  10"’ 


Time,  sec 
(a) 
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(c) 

Figure  2.4  Central  responses  of  the  plate  at  (a)  (0=3  rad/s,  (b)  co=100  rad/s,  (c)  co=200  rad/s 


2.4  Summary 

A  Co-Rotational  formulation  for  static  and  dynamic  analysis  of  geometrically  non-linear  plates  is 
presented  using  DKT  triangular  elements.  A  number  of  numerical  examples  are  provided  and  the 
results  are  compared  with  other  methods  of  computation.  The  computation  in  each  loading  was 
obtained  by  using  a  single  load  increment  i.e.  the  load  for  which  the  solution  is  sought  is  applied  in 
full  at  a  time  and  the  iteration  was  performed  within  this  load  until  the  difference  of  the  applied  full 
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load  and  the  resisting  forces  becomes  minimum.  The  computation  can  also  be  performed  by  using  as 
many  load  steps  as  desired  without  using  the  full  load  at  a  time.  But  the  number  of  iterations  was 
found  to  be  much  less  for  the  former  case.  For  proper  convergence  at  higher  load  levels  a  relaxation 
factor  has  been  introduced.  For  static  analysis  this  fomulation  produces  the  closest  results  of  all  the 
methods  to  the  exact  solution.  There  are  some  discrepencies  in  responses  between  the  present  method 
and  other  methods  in  dynamics  analysis  of  plates  which  needs  to  be  investigated  in  future. 


2.5  Appendix 


2.5.1  Determination  of  Rotational  Transformation  Matrix 

1.  Determine  the  local  co-ordiantes  of  the  imdeformed  triangular  element  ),  (Xj.Fj ) 

and  Xj ,  Yj  from  the  input  global  co-ordiantes  and  determine  the  co-ordinates  of  Q 


= 


X,  +X3  VFY3 


(A2.1) 


and  the  angle 


p  =  tan"' 


'Mi' 


(A2.2) 


2.  Define  local  co-ordinate  system  by  identifying  the  X'-Y‘‘  plane  by  passing  through  all  the  three 
displaced  nodes  1 , 2,  and  3 .  Determine  the  co-ordinates  of  Q  as 


x**  v**  7^  = 


'"x^-i-x^  y2+y|  4+4^ 


(A2.3) 


2  2  2 

/ 

3.  The  undeformed  element  is  so  placed  on  the  displaced  plane  that  the  first  node  of  the 
undeformed  element  coincides  with  the  first  node  of  the  deformed  element  (origin  of  5^-Y  -Z  co¬ 
ordinate  system),  i.e  (xf,Y,'*)=  0,0  .  (x2,Y2)and  Xj.Yj^  can  be  obtained  by  coordinate 
transformation  Establish  the  vector 


V, 


Qy 


V,.q.=Vq.  =  X-  -  X?  y“  -y?  -z?  = 

Similarly  vector Vj, 2'  =  '^2  y\  ~  y?  4~4} 

A  vector  V2d  in  the  direction  Z"*  is  given  by  the  cross  product  =  Vi.q.  X  Vj.2. 


V, 


Qz 
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yj-y?  4-4  -  yl-y?  4-4 

II 

(z J  - zf  Ixj -xf  )-(xJ  - xf  Izj  - z?) 

(xj  -  X?  lyt  -  yf )-  (x^  -  xf  )(yj  -  yf ) 

Or,  Vgi  -  V2X 

Vzy  V2, 

4.  Now  the  vector  perpendicular  to  Z**  - 1‘  -  Q'  plane  is  given  by  Vj^d  =  V^d  x 

Vj^d  =  V^y  Vq2  ~  V^z  Vqy  “  V2xVqz  ^Zx^Qy  “  ^Zy  Vqx 

5.  Determine  the  magnitudes  of  the  vectors  Vpd  jV^d  and  Vj^d  as 

Rqd  “  ^|^Qx  "^^Qy  "*"^Qz 

R-z-  =  -y/vlx  +V|y  +^Zz 

RR-=VVfo+V|y+V|,  (A2.4) 

6.  Let  (l,,nii,n,)  denotes  the  direction  cosines  of  the  vectors  Vqj  .  Similarly 


IjjUljjnjandlj.mjjnj  are  the  direction  cosinses  of  the  vectors  V^d  and  respectively.  The 
direction  cosines  are  given  by 

Vn 


1  n,  - 

-  nil  - 

Kgd 


1,= 


V, 


Rx 


Qy  „  _ 

- ,  Ui  — 

Rq<j  Rqd 

Vp 


R 


= 


,  n^  =■ 
R,d  R, 


U  = 


R^d 


Zx  ^  _Vzy  ^  _V,, 

.  1113  — »  “3  -■ 


R. 


R. 


(A2.5) 


The  direction  cosines  of  the  new  co-ordinate  system  are  given  in  matrix  form  as 


■q"' 

1. 

mi 

ni  ' 

x'' 

R** 

= 

12 

m2 

“2 

y' 

Z'* 

^3 

m3 

z' 

(A2.6) 


7.  Now  rotate  the  coordinate  system  around  Z*  by  the  angle  P,  and  determine  die  transformation 


matrix  T,  from 
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X'*' 

COsP 

sin  p 

o' 

pi 

mi 

ni  ' 

x' 

yd 

= 

-sin  P 

cosP 

0 

^2 

m2 

n2 

y' 

z** 

0 

0 

1 

;3 

m3 

ns. 

z' 

■V- 

T. 


(A2.7) 


The  elements  of  the  transformatin  matrix  is  given  by  aij  as 

[tJ  = 


a.i 

^12 

^13 

^2! 

a  22 

^23 

.^31 

^32 

^33  _ 

(A2.8) 


8.  There  are  five  degrees  of  freedom  at  each  node  of  the  element;  three  translational  components 
and  two  rotational  components.  Hence  the  complete  transformation  matrix  of  the  element  would  be 


[t]= 


0 


where  [t^]  = 


“21 

a3i 

0 

0 


“12 

a22 

a32 

0 

0 


0  X 


(A2.9) 


cj 


a, 3  0 

^23  0 


*33 


0 

0 


*21 


12 


*22  J 


Steps  in  getting  the  transformation  matrix  T  : 

15x15 

1.  Determine  the  local  coordinates  of  Q  from  equation  (A2. 1); 

2.  Compute  the  angle  P  (location  of  the  median  of  the  element  fi'om  the  local  x-axis)  for  each 
element  from  equation  (A2.2); 

3.  Compute  the  directions  cosines  of  the  new  coordinate  systems  fi'om  equation  (A2.5); 

4.  Compute  the  elements  of  the  trasformation  matrix  for  three  translational  components  of 
displacement  from  equation  (A2.7); 

5.  Determine  the  complete  15x15  transformation  matrix  T  form  equation  (A2.9). 


2.5.2  Derivation  of  Stiffness  Matrix  for  Plane  Stress  Condition 

Element  displacement  is  represented  in  terms  of  nodal  displacement  as 


55 


Computational  Mechanics  Laboratory,  University  of  Arkansas  at  Fayetteville  Selvam,  Roy,  and  Qu 

u®  =  a,  +a2X  +  0C3y  or,  u®  =  1  x  y  aj  (Xi  013  ^  (A2.10) 

Similary, 

v‘=Pi  +p2X  +  p3y  or,  V®  =  1  X  y  p,  p2  p3  ^  (A2.11) 

Local  coordinates  of  nodes  1,  2  and  3  are  (xi,yi),  (X2,y2)  and  (X3,y3),  respectively.  From  equation 
(A2.1),  we  have 

,  =a,  H-ttjXi  H-ttj  , 

2  =  a,  +a2X2  +063  2 

3=a,+a2X3+a3  3 

which  can  be  written  in  matrix  form  as 

1 

(A2.12) 


_ 1 

1 - 

1 _ 

- 1 

P 

_ 1 

^2 

= 

1  X2  y2 

“2 

"3. 

.1  X3  y3 

or, 


a, 

1 

_ 1 

-1 

U| 

“2 
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1  X2  y2 

U2 

I 

I 

>> 

_ 1 

The  coefficient  matrix  can  be  obtained  as 


a, 

1 

_ 1 

“2 

II 

U2 

1 

_ 1 

.«3. 

(A2.13) 


or 


where. 


a. 

“2 

1 

"lAl 

“3. 

X3y,-x,y3  x,y2-X2y, 


y2-y3 


X3  -X2 


y3-yi 


X,-X3 


yi-ya 


X2-X, 


Ul 

U2 

.'^3. 

|A|  = 


1  X, 
1  X2 
1  X, 


y2 

y3 


=(x2y3  -X3y2)-(x,y3  -X3y,)+(x,y2  -X2y,) 

=2x  area  of  the  triangle=2A 
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Hence  equation  (A2.10)  can  be  expressed  in  terms  of  nodal  co-ordinates  and  nodal  displacement 
components  as 


u'=[l  X  y]j^ 


|A| 


XjYi-XiYj  XiYz-Xjy, 


y2-y3 


X3-X2 


ya-yi 


X,-X3 


yi-y2 

X,  -X, 


1 

c 

_ 1 

«2 

- 1 

«n 

_ 1 

u®  =  N,  x,y  N2  x,y  N3  x,y  Uj  U2  U3 
or,  u'  =  N  x,y  [u] 


(A2.14) 


N,(x,y)  N2(x,y)  N3(x,y)  0  0  0 

0  0  0  N,(x,y)  N2(x,y)  N3(x,y)J 


u, 


u. 


u. 


LV3J 


N,{x,y)  0  N2(x,y)  0  N3(x,y)  0 

0  N,(x,y)  0  N2(x,y)  0  N3(x,y) 


u, 


u. 


LV3J 


(A2.15) 


ep  = 


L^xyj 


a 

dx 


^  0 


0 


ay 

i  1 

[ay  axj 


Strain-displacement  relationship; 


Ep  =  Bp  ap 


(A2.16) 


where 


KI4 


Yi-Vi  0  y3-yi  0 


0  X3  -  X2 


0 


yi  -y2 

0 


X3-X2  y2-y3  x,-X3  y3-y,  X2-X,  y,-y2 


and 
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Kp  =  j  B/DBptdA 

A 

=  [B,fD[B,]tA 

For  plane  stress  condition,  the  elasticity  matrix  is  given  by 


1 

V 

0 

D  = 

ZZ  CN 

1 _ 

du 

^22 

o' 

0 

\> 

1 

0 

0 

0 

d33. 

0 

0 

l-\) 

2  J 

(A2.17) 


(A2.18) 


Note:  The  elements  of  the  stiffiiess  matrix  for  plain  strain  condition  are  given  in  Reference  2,  pp.69. 


2.5.3  Derivation  of  Stiffness  Matrix  for  Piate  Bending 

In  choosing  a  function  for  the  lateral  displacement  w,  the  polynomial  must  include  nine  unknown 
constants  e.g. 

w  =  a,  +  ttjx + ajy  +  a4X^  +  ccj  y  V  a^xy + ayx^  +  agX^y +095^ 

Comparing  this  to  a  full  cubic  polynomial,  it  may  noted  that  one  term  has  been  omitted,  namely  xy 
term.  Many  researchers  have  tried  several  methods  to  overcome  this  difficulty.  For  example  the  ten 
terms  of  the  full  cubic  expression  can  be  retained  in  the  expression  for  w,  two  of  the  coefficients  being 
specified  to  be  equal,  e.g.  ttg  =  09.  However,  if  this  done,  the  [A]  matrix  in  the  expression  A2.4 
becomes  singular  for  certain  orientation  of  triangular  element,  e.g.  when  two  sides  of  the  isosceles 
triangles  are  parallel  to  x  and  y-axes. 

Derivation  of  element  stiffness  matrix  using  area  coordinates 
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a,  +b,x+c,y 

L  =— ^ - - - and  similarly 

1  2A 


aj+bjX  +  Cjy  a3+b3X  +  C3y 

^2  -  2 A  2 A 


(A2.22) 


where 


^2  =  ^3X1  -  YaX,  .bz  =  Ys  -  Yi ><^2  =  Xi  -Xj 
aj  =x,y2-YiX2.b3  =  y,-Y2.C3  =X2-x, 

Total  displacement  in  the  elemenh=rigid  body  movement  (no  curvature)+simply  supported  case 
(No  nodal  translation) 

W  =  w*+w®*  (^-23) 

Since  no  curvatures  as  set  up  in  the  case  of  rigid  body  movement  case,  the  lateral  displacement  must 
be  linear  function  of  x  and  y 

=  w  jLj  +  W2L2  +W3L3  (A2.24) 

Consideraing  simply  supported  case,  the  element  has  only  two  rotaional  degrees  of  freedom  at  each 
node,  i.e.  and  0y.  The  element  displacment  vector  may  be  written  as: 

[D“]=e“,  e“  9S  9“  9S  6“^  (A2.25) 


The  corresponding  force  vector  F®®  contains  six  terms  and  the  stiffness  matrix  K®®  is  then  a  6X6 
matrix. 

ay  X 


(A2.26) 


Substituting  for  w“  from  equation  (A2.23) 

9..=e,+^Md9';=9,-^ 

Substituting  the  expression  assumed  for  w*  in  equation  (A2.24)  into  equation  (A2.26),  the  slopes  at  a 
point  in  the  simply  supported  case  can  be  expressed  in  terms  of  total  slopes  at  that  point  and  the  lateral 

displacement  as  follows: 

Gx  =0*  +— (w,c, +  W2C2+W3C3) 


0y  =®y  -^(w,b,+W2b2+W3b3) 

y  y  ^ 


(A2.27) 
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The  slopes  can  be  determined  at  each  node.  For  example  at  node  1, 

0x1  =0xl+^(w,C,+W2C2+W3C3) 

A 

0yl  =0yl  -^(W,b,+W3b2+W3b3) 

Similar  expression  can  be  obtained  for  other  two  nodes  and  the  complete  relationship  between  the 
nodal  rotations  of  the  simply  supported  element  and  those  of  the  complete  element  may  be  written  in 
matrix  form  as 


[d®“]= 


! 

<x> 

1 

cl 

2A 

0 

^2 

0 

0 

^3 

0 

0  ] 

CD 

-bl 

0 

2A 

-b2 

0 

0 

-b3 

0 

0 

/\SS 

^x2 

1 

cl 

0 

0 

c2 

2A 

0 

c3 

0 

0 

^y2 

~2A 

-bl 

0 

0 

-b2 

0 

2A 

-b3 

0 

0 

qss 

cl 

0 

0 

c2 

0 

0 

c3 

2A 

0 

QSS 

-bl 

0 

0 

-b2 

0 

0 

-b3 

0 

2A  J 

This  can  be  summarized  as: 

D®®  =[t]  D® 

From  the  princple  of  virtual  displacement 
pess  pess  _  j^e  pe 

From  equations  (A2.29)  and  (A2.30), 

F®  =  [t]^  F®*® 

Again, 

F®ss  _  j^ess 

From  equations  (A2.32)  and  (A2.29),  it  can  be  shown  that 
Kg  =  K'  =  T  ^  K'**  T 


w, 


xl 


Uy, 

W2 

0x2 

Wj 

0x3 


(A2.28) 


(A2.29) 

(A2.30) 

(A2.31) 

(A2.32) 

(A2.33) 


Derivation  of  Stiffness  Matrix  for  Simply  Supported  Element 
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Lateral  displacement  at  a  point  on  the  simply  supported  element  in  terms  of  the  nodal 
displacement  is  chosen  as  follows: 

W“  =N,,e“  +N,,e“  +N^e5,  +N,,ej  (a2.34) 

Each  individual  Nx  and  Ny  term  in  the  expression  represents  the  sh^e  fimction.  The  slopes  0®®  and 

oSS 

0®®  at  each  point  are  obtained  from0j^  =  ~z  and  0^  —  — • 


dx 


QSS  qSS  ,  ^^05S  qSS  ^ 

-  ay  ay  ay 


ILqss  .  rvSS  .  /xSS  .  ^^x3  qss 


N 


-0„.  +■ 


y30ss 


3y  ”y2  ^  3y  ^*3  ^  3y  y3 


0“  =- 
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’3N 


xl  0SS  _J_ 


aN 


ax  ''y3 


aN 


aN 


-0®®  + 


0"l+- 


aN. 


aN 


y2  nss  .  ‘^^''x3  r^ss  .  "‘'y3  yyss 


-0’*l  + 


-0® 


The  individual  shape  functions  must  be  chosen  so  that  they  satisfy  boundary  conditions  at  the  nodes  of 
the  simply  supported  element  and  constant  strain  criterion. 


f  0  1 

Nxl==b3  L1L2  +  — L1L2L3 


Nyl  —  C3  L1L2  +  L1L2L3 


h( 

—  C2^LiL3  +  — 


2  1 

LjLj  "H  “  L1L2L3 


L1L2L3 


(A2.35) 


The  expression  for  other  shape  functions  can  be  obtained  by  changing  the  subscripts  in  cyclic  order. 


1  \ 

=  bj|  L2L3  “^”’^L|Lf2L3 


L2L3  +  TL1L2L3  1“C3  L2L1  +';:7LiL2L3 


Nx2=b,|^ 

Ny2=C.|^ 

Nx3=b2j^LlL,+i 

Ny2  =C2^L3Li  +  — 


^  2  1 

L2L1  "^~L|L2L3 


) 

) 


2  ^ 

L|Lf2L3  — bj  L3L 2  I^2^3 


LIL2L3  L3L2 


iL,L2L3j 


The  curvatures  set  up  in  the  simply  supported  element  may  now  be  related  to  nodal  displacements  as 
follows; 

,  f  :32„,ss  -'2  "  -'2 

[eix,y)]= 


a^w 


ax^ 


a'w“ 


ay^ 


axay 
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SS  "dSS  tvSS 

8  x,y  =  B  D 

where  the  matrix  B*®  3X6  matrix  and  can  readily  be  determined  by  substituting  the  expression  for  w'- 
from  equation  (A2.34). 


[b“]=  - 


aX2 

a^Ny, 

a^N,3 

a^N,3 

ax' 

ax' 

ax' 

ax' 

ax' 

ax' 

ax. 

a^N,. 

a^Nx2 

aX2 

a^Nx3 

a^N,3 

ay' 

ay- 

ay- 

ay- 

ay- 

ay- 

,aX. 

2a'N„ 

ja-N, 

ja-N., 

axay 

axay 

axay 

dxdy 

axay 

dxdy 

(A2.36) 


Finally  the  stiffhcss  matrix  for  the  simply  supported  element  can  be  obtained  from 

K*®  =|J  B®®  B®®  tdxdy  (A2.37) 

where  K®®  is  a  6  x  6  matrix.  Having  thus  found  the  stiffness  matrix  for  the  simply  supported 

element,  the  stiffiiess  for  the  actual  element  K®  can  be  obtained  from  equation  (A2.33). 

Steps  in  deriving  the  element  stiffiiess  matrix  for  plate  bending: 

1.  Determine  the  local  coordinates  of  elements  using  the  transformation  matrix  as  described  in 

appendix  A-1; 

2.  Compute  the  values  of  aj,  bi  and  Cj  from  equation  (A2.2 1)  and  determine  the  area  coordinate  values 
Li,L2  andLs; 

3.  Compute  the  values  of  shape  functions  Nxi,  Nyi  etc.  from  equation  (A2.35); 

4.  Establish  the  B”  matrix  in  equation  (A2.36); 

5.  Evaluate  the  product  B®®^EB®®  ; 

6.  Integrate  this  product  as  in  equation  (A2.37)  to  get  the  stiffness  matrix  for  the  simply  supported 
element 

7.  Evaluate  T  matrix  as  given  in  equation  (A2.28); 

8.  Evaluate  the  required  element  stiffiiess  matrix  Kg  from  equation  (A2.33); 

In  this  computation  a  closed  form  DKT  element^  is  used  for  bending. 
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MOVING  GRID  AND  DEFORMING  MESH 


3.1  Introduction 

3.2  Description  and  Evaluation  of  Different  Algorithms 

3.3  Theory  of  the  Deforming  Mesh 

3.4  Numerical  Simulations 

3.5  Siunmary 


3.1  Introduction 

As  we  know,  fluid  equations  are  derived  using  the  Eulerian  motion  description  (stationary  grid 
with  respect  to  time)  and  the  stmctural  equations  derived  using  the  Lagrangian  description  (grid  moves 
with  the  displacement  of  the  physical  medium).  To  accurately  model  problems  with  moving  and 
deforming  boundaries,  such  as  in  aeroelastic  analysis,  requires  that  the  fluid  grid  conform  to  the 
changing  shape  of  the  boundary.  To  accommodate  the  movement  of  the  fluid  grid  in  an  aeroelastic 
code  either  a  new  grid  must  be  generated  at  each  time  step  (dynamic  regridding)  or  the  differential 
equations  for  the  fluid  must  be  modified  to  allow  the  movement  of  the  existing  fluid  grid  (ALE, 
Corotational,  Space-Time).  Once  a  method  is  found  for  modifying  the  fluid  equations  to  allow  for  grid 
movement  a  method  must  be  chosen  to  move  or  update  the  fluid  grid  to  the  new  position  of  the 
structure. 

Grids  are  usually  moved  or  updated  with  the  same  constraints  used  in  the  programs  used  to  create 
the  initial  grid.  Such  constraints  control  the  size  of  the  grid  cells,  their  distribution,  and  the 
orthogonality  of  the  cell  intersections  and  boundaries.  The  goal  of  any  grid  updating  or  moving 
method  is  to  maintain  a  desirable  grid  throughout  the  duration  of  the  program.  If  at  any  point  the  grid 
fails  to  meet  the  required  level  of  quality  a  new  grid  must  be  generated  that  meets  the  desired  quality. 
Generating  a  new  grid  or  regridding  requires  the  incorporation  of  a  grid  generation  program  into  the 
aeroelastic  solver  and  the  interpolation  of  the  values  from  the  old  grid  to  the  new  grid.  This  process 
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takes  time  and  can  introduce  errors  due  to  the  interpolation  process.  However,  such  a  program  is  often 
needed  for  aeroelastic  problems  with  large  deflections  and  deformation. 

Ideally,  a  grid  would  be  moved  in  such  a  way  that  would  maintain  grid  quality  under  large 
deflection  and  deformation,  introduce  minimal  errors,  and  be  computationally  efficient.  Just  like  any 
other  numerical  method  you  have  to  give  in  one  of  these  areas  to  get  in  another.  One  way  to  increase 
computational  efficiency  was  addressed  by  Donea'  in  that  since  the  computational  domain  or  the 
boundary  of  the  fluid  grid  usually  extends  relatively  far  away  from  the  structure  and  its  movement, 
only  the  portions  of  the  fluid  grid  near  the  structure,  are  required  to  be  moved  to  accommodate 
deflection  and  deformation.  The  less  grid  points  that  are  moved  at  each  time  step,  the  more  efficient 
the  code  will  be^  However,  grid  quality  may  suffer  if  the  area  allowed  to  move  is  not  large  enough  to 
allow  a  smooth  transition  from  the  fixed  points  and  the  displaced  structure. 

In  this  chapter  the  different  methods  used  to  modify  the  fluid  equations  to  allow  the  movement  of 
a  structure  within  the  fluid  domain  are  described.  For  the  ALE  method,  various  methods  used  to 
update  or  move  the  fluid  grid  to  conform  to  a  changing  structural  boundary  are  also  described  and 
evaluated.  The  details  of  the  implementation  of  the  Trans-Finite  Interpolation  scheme  for  the 
deforming  mesh  is  provided.  Numerical  simulation  is  also  included  to  show  the  features  of  this 
scheme. 


3.2  Description  and  Evaluation  of  Different  Algorithms 

3.2.1  Dynamic  Regridding 

The  dynamic  regridding  method  is  implemented  by  simply  forming  a  new  grid  at  each  time  step  to 
conform  to  the  shape  of  the  stmcturel  The  method  can  be  summarized  by: 

•  Create  an  initial  fluid  grid. 

•  Calculate  the  fluid  pressure  on  the  stmcture. 

•  Calculate  the  stmctural  response. 

•  Generate  a  new  grid  around  the  new  position  of  the  structure. 

•  Obtain  the  values  of  the  flow  variables  for  each  grid  point  of  the  new  grid  by  interpolating  the 
values  of  the  flow  variables  from  nearby  grid  points  of  the  old  grid. 

The  new  grid  needed  at  each  time  step  can  either  be  generated  manually  (i.e.  a  person  with  a  grid 
generation  program)  or  by  an  automatic  grid  generation  program  incorporated  into  the  aeroelastic 
program.  A  simple  aeroelastic  analysis  can  require  several  thousand  time  steps.  This  makes  it 
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impossible  to  generate  the  grid  manually.  Programs  with  the  grids  generated  automatically  require  a 
very  robust  generation  program  to  ensure  grid  quality  since  the  grids  generated  automatically  cannot 
be  viewed  for  quality  and  errors  during  the  execution  of  the  program. 

Compared  to  newer  methods  for  moving  boundaries,  dynamic  regridding  is  inefficient  due  to  the 
amount  of  time  needed  to  generate  new  grids.  The  method  also  introduces  greater  errors  due  to  the 
interpolation  of  data  from  the  old  grid  to  the  new. 

3.2.2  Corotational  method 

Among  the  different  methods  developed  to  allow  solvers  to  handle  moving  fluid  grids  rs  the 
corotational  method.  Farhat  and  Lin'*  described  the  equation  of  motion  of  the  fluid  in  moving  frames  of 
reference,  which  are  attached  to  specific  nodes  of  the  structure.  The  result  is  an  implicit  generation  of  a 
structure  attached  corotational  fluid  grid.  Using  this  scheme  the  grid  does  not  need  to  be  updated  for 
rigid  body  motion.  Instead  of  recalculating  the  new  position  of  the  grid  points,  the  entire  coordinate 
axis  system  is  moved  or  rotated  in  physical  space  as  seen  in  Figure  3.1. 


Figure  3.1  Moving  frame  of  reference 


The  result  of  this  procedure  is  that  the  coordinates  of  the  individual  grid  points  do  not  have  to  be 
recalculated  at  every  time  step.  The  grid  is  moved  in  the  physical  domain  but  remains  stationary  in  the 
computational  domain^  The  fluid  equations  are  modified  by  adding  terms  that  describe  the  grid  and 
axis  movement  on  the  right  hand  side  of  the  equations.  This  method  was  originally  designed  for  rigid 
body  motion,  but  it  can  be  modified  to  account  for  small  structural  deformations.  Lin  described  ways 
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to  implement  the  corotational  approach  to  flexible  configurations  by  implementing  different  Jacobian 
updating  schemes. 

The  attraction  of  this  method  is  its  ability  to  move  the  entire  coordinate  system  with  the  movement 
of  the  stmctuie.  No  updating  is  needed,  so  the  same  grid  can  be  used  for  all  time  steps.  The  downside 
of  this  method  is  the  difficulty  in  applying  the  unusually  large  amount  of  terms  added  to  the  fluid 
equations  to  account  for  movement.  If  small  structural  deflections  are  to  be  accounted  for  even  more 
terms  are  added. 

3.2.3  Arbitrary  Lagrangian  Eulerian  method 

A  method  presented  by  Hirt,  Amsden,  and  Cook®,  called  the  Arbitrary  Lagrangian  Eulerian 
method  or  ALE  for  short,  suggests  a  way  of  modifying  the  fluid  equation  in  a  way  that  would  reflect 
Eulerian,  Lagrangian  or  an  arbitrary  combination  of  the  two  descriptions  of  motion.  This  method, 
originally  developed  for  finite  difference  equations,  was  later  modified  for  use  with  finite  elements' . 

The  ALE  method  takes  advantage  of  the  good  properties  in  the  Lagrangian  and  Eulerian  systems 
by  allowing  a  specified  region  of  the  fluid  grid  to  move  to  conform  to  structural  displacements.  Yet 
allow  the  outer  regions  of  the  grid  to  remain  motionless.  Figure  3.2  shows  the  location  of  the  ALE 


Pure  Eulerian  Coordinates 

Arbitrary  Lagrangian 
Eulerian  Coordinates 


Pure  Lagrangian  coordinates 


Figure  3.2  Body  conforming  grid  with  ALE  region 

portion  of  the  fluid  grid.  In  this  figure  the  different  descriptions  of  motion  that  correspond  to  a  specific 
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grid  region  are  shown.  The  blue  shows  the  Eulerian  coordinates,  the  green  is  Lagrangian,  and  the  red 
is  the  ALE. 

In  this  method  the  equations  for  the  structure  and  the  outer  fluid  region  remain  unchanged  with 
T  agrangian  and  Eulerian  descriptions  respectively.  For  the  ALE  portion  of  the  fluid  grid  the 
convective  flux  terms  of  the  fluid  equations  are  modified  to  account  for  the  change  in  velocity  of  the 
moving  grid.  This  is  done  by  taking  into  account  the  fluid  grid  velocity  and  subtracting  it  firom  the 
calculated  fluid  velocity  at  each  node.  The  equation  reprinted  below  is  an  example  of  the  ALE  form  of 
the  momentum  equation.  The  highlighted  portion  of  the  equation  shows  the  grid  velocity  Vj  being 


+1 


Wu+'JlX 


subtracted  from  the  fluid  velocity  Uj.  The  grid  velocity  is  calculated  by  taking  the  change  in  position 
of  the  grid  with  respect  to  time  from  one  time  step  to  the  next.  This  velocity  will  vary  depending  on 
the  location  within  the  ALE  grid.  At  the  surface  of  the  structure  the  grid  velocity  will  be  equal  to  the 
velocity  of  the  structure  and  decrease  with  distance  from  the  structure  until  it  is  equal  to  zero  at  the 
outer  edge  of  the  ALE  grid  as  seen  in  Figure  3.3.  The  result  is  a  modified  fluid  equation  that  represents 


Figure  3.3  Grid  velocities  for  ALE  method 

pure  l  agrangian  motion  at  the  structure  surface,  pure  Eulerian  motion  description  at  the  outer  edge  of 
the  ALE  grid,  and  a  mk  between  Lagrangian  and  Eulerian  motions  descriptions  in  between. 
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The  ALE  method  is  one  of  the  most  widely  used  methods  for  problems  involving  moving 
structures.  One  advantage  to  this  method  is  the  ability  to  specify  how  much  of  the  fluid  grid  is  allowed 
to  move.  Limiting  the  size  of  the  ALE  region,  limits  the  amount  of  grid  that  needs  to  be  updated  to 
correspond  to  the  new  position.  This  saves  computational  time  and  assures  grid  quality  for  the 
unmoved  portion  of  the  grid.  The  method  can  handle  rigid  body  problems  as  well  as  problems 
involving  structural  deflection.  Unlike  the  corotational  mefliod,  the  ALE  method  is  also  suitable  for 
problems  where  the  outer  perimeter  of  the  fluid  grid  is  not  allowed  to  move  due  to  a  boundary 

condition.  However  the  method  does  have  a  disadvantage. 

The  disadvantage  of  the  ALE  method  is  that  it  requires  that  the  fluid  grid  be  updated  after  each 
time  step  to  correspond  to  the  new  position  of  a  structure.  The  quality  and  efficiency  of  the  ALE 
approach  are  dependent  on  the  grid  updating  procedure  used.  If  an  updating  procedure  used  fails  to 
maintain  grid  quality,  then  the  grid  must  be  regenerated,  which  is  computationally  expensive. 
Therefore  a  grid  updating  procedure  should  be  chosen  carefully.  Several  grid  updating  procedures 
include  the  Trans-Finite  Interpolation  (TFI),  a  Cubic  Blending,  a  dynamic  mesh,  and  a  rigid  grid 
method  are  discussed  below. 

3.2.3.1  Trans-FInlte  Interpolation  (TFI) 

A  commonly  used  method  for  grid  generation,  called  Trans-Finite  Interpolation  (TFI),  was  used 
by  Guruswamy^.  In  this  paper  Guruswamy  developed  an  algebraic  method  using  TFI  for  aeroelastic 
configuration  adaptive  grids.  TFI  controls  grid  quality  by  connecting  grid  lines  perpendicular  to  the 
structure  and  distributing  the  nodes  exponentially  away  from  the  structure  to  the  outer  boundary.  The 
outer  boundary  is  also  located  a  significant  distance  away  from  the  moving  stmcture.  The  outer 
portions  of  the  grid  are  restrained  from  moving  by  shearing  the  grid  along  the  outer  boundary.  Grid 
lines  are  prevented  from  overlapping  by  preserving  the  arc-length  distribution  between  nodes,  and  by 
using  a  simple  linear  distribution  of  translational  displacement. 

At  each  time  step  the  deformed  shape  of  the  structure  is  calculated  by  the  displacement  vector 
{d}=[0]{q},  where  [O]  is  the  modal  matrix  and  {q}  is  the  generalized  displacement  vector.  The  new 
positions  of  the  grid  points  are  calculated  using  the  positions  of  the  grid  points  from  previous  time 
steps.  This  is  achieved  using  a  first  order  backward  difference  scheme.  The  grid  points  are  distributed 
along  the  grid  lines  in  the  radial  direction  using  a  spacing  function.  An  example  of  this  updating 
scheme  can  be  seen  in  Figure  3.4. 
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Initial  grid  Displaced  grid 

Figure  3.4  An  example  of  the  TFI  updating  scheme* 


Shankar  and  Ide’,  Gaitonde  and  Fiddes*,  and  Gaitonde,  Jones,  and  Fiddes*®  used  a  general 
structured  grid  method  based  on  TFI  for  rapid  grid  generation  at  each  time  step.  Selvam,  Visbal,  and 
Morton"  used  TFI  method  to  update  grid  points  in  an  aeroelastic  solver  for  non-linear  panel  flutter. 
The  advantages  of  the  TFI  method  is  that  it  is  efficient  and  relatively  simple  to  apply.  A  disadvantage 
is  that  mesh  orthogonality  at  deforming  surfaces  cannot  be  guaranteed  with  TFI  dining  moderate 
deflections".  Orthogonality  at  the  surface  is  desired  for  simulations  with  a  high  Reynolds  number  and 
viscous  flow.  In  some  cases  TFI  can  limit  solution  accuracy  and  stability. 


3.2.3.2  Cubic  Blending 

Another  algebraic  method  was  developed  by  Melville,  Morton,  and  Rizzetta".  This  method  is 
similar  to  TFI  in  that  it  is  based  on  redefining  the  grid  lines  normal  to  the  surface  but  it  assures  grid 
orthogonality  near  a  deforming  surface  for  deflections  and  rotations.  This  method  also  allows  the 
option  of  fixing  the  outer  region  so  the  grid  overlap  connectivity  remains  unchanged.  The  name  cubic 
blending  comes  from  the  procedure  used  to  combine  an  old  grid  and  a  calculated  reference  grid 
together  to  form  a  new  grid. 

The  procedure  for  this  method  starts  by  solving  the  fluid  equations  and  transferring  the  pressure  to 
the  structure.  Then  after  the  structural  equations  have  been  solved  and  the  structure  is  moved,  the  new 
grid  is  calculated  by  first  calculating  the  transition  and  rotation  of  each  surface  node  by  the  equations: 

Translational  displacements:  Ax  =  x,. ,  -  x, ,  Ay  =  y,. ,  -  y,_, 

„  s  s  sxs  -e^ 

Rotational  displacements:  cosa,  —  pjp|  smp 

Where:  s  is  the  original  surface  vector  from  i-1  to  i+1 
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5'  is  the  displaced  surface  vector 

Then,  according  to  the  displacement  of  the  surface  node,  the  normal  grid  line  is  moved  as  a  rigid  body 
to  form  a  displaced  grid  line,  which  serves  as  a  reference  grid  line.  These  referenced  grid  lines  are 
defined  by  the  equations; 

-^u)cos0,  -(y,.,. -y,.,,)sin0, 

yUr  =  yi.i  +  ^y,- + 1.  j  ~ 

Blending  the  old  grid  line  and  the  reference  grid  line  together  forms  the  new  grid  line.  Morton  suggest 
blending  the  grid  lines  in  arc-length  space  rather  than  in  computational  space  by  using  the  arc-length 
equation: 

^ -J {xjjc  ~ ^i.it-1 )  '^{yijc'^ yi,k-\ ) 

k=2 

Where  the  arc-length  of  the  starting  node  is  5, ,  =0.  Cubic  blending  with  a  zero  slope  at  the  end 

points  ensures  wall  orthogonality  is  maintained,  and  that  there  is  a  smooth  gnd  translation  in  die  outer 
regions  of  the  grid.  The  blending  equation  is: 


The  variable  j  max  is  the  last  node  along  each  normal  line  that  is  allowed  to  deflect.  The  new  grid  can 
now  be  formed  by  ^plying  the  blending  function  to  the  reference  grid  and  the  old  grid  by  the 
equations: 


Figure  3.5  Displaced  grid  on  a  circular  cylinder  using  the  cubic  blending  method'^ 
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The  advantage  of  this  method  is  that,  unlike  TFI,  cubic  blending  allows  the  grid  lines  to  intersect 
perpendicular  to  the  structure  and  the  outer  boundary.  By  doing  this  the  quahty  of  the  grid  can  be 
maintained  under  larger  deflections.  The  method  is  also  applied  in  a  relatively  simple  manner.  One 
disadvantage  is  that  cubic  blending  may  require  more  grid  points  to  be  allocated  to  move  than  TFI. 
This  is  done  to  ensxue  grid  quality. 

3.2.3.3  Dynamic  mesh 

A  unique  way  of  updating  the  fluid  grid  was  developed  by  Batina"  for  an  unstructured  mesh.  This 
method  uses  the  ALE  method  in  that  it  allows  the  outer  regions  to  be  held  fixed  and  a  region  to  be 
specified  where  the  fluid  grid  is  allowed  to  move  and  deform  to  conform  to  the  deflected  or  deformed 
structure.  In  this  method  Batina  models  each  edge  of  a  triangular  element  in  the  fluid  grid  by  a  sprmg, 
where  the  spring  stiflhess  is  the  inverse  of  the  length  of  the  edge  of  the  element  and  calculated  by  the 

equation  kr.  =  \l  ^{y ,-yif  ^  Figure  3.6  shows  an  example  of  a  triangular  element 

modeled  as  a  spring.  An  example  of  a  network  of  triangular  springs  can  be  seen  in  figure  3.7.  In  this 
figure  the  blue  indicates  the  stationary  Eulerian  grid  points,  and  the  red  the  ALE  grid  points. 


Figure  3.6  Element  modeled  as  a  spring  system 
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Figure  3.7  Illustration  of  a  fictitious  spring  system  about  an  airfoil. 


The  procedure  for  this  method  starts  by  solving  the  fluid  equations  and  transferring  the  pressme  to 
the  structure.  Then,  the  structural  displacement  is  calculated  and  the  structure  is  moved.  The  new 
mesh  is  found  using  a  predictor-corrector  procedure.  At  each  time  step  when  the  structure  deflects,  the 
grid  points  on  the  outer  boundaries  are  held  fixed  and  the  grid  points  on  the  body  are  prescribed  by  the 
structural  deflection.  Moving  the  nodes  along  the  structure  to  their  new  location  causes  tension  forces 
in  the  network  of  the  springs.  Summing  these  forces  at  each  point  in  the  x  and  y  direction,  gives  static 
equilibrium  equations.  These  equations  are  then  solved  for  the  nodal  displacements,  5,,  and  .  The 
nodal  displacements  are  predicted  by  a  linear  extrapolation  of  displacements  from  two  previous  time 
steps. 

Nodal  displacements:  S^=25",-S"~'  =2d”, -S”,' 

These  displacements  are  then  corrected  using  Jacobi  iterations  of  the  static  equilibrixun  equations. 

\  k  6  ,  y' 

sr/i+l  _  ^n+\  _  ^  ym 

■  s*- 

Then  summing  for  all  the  edges  of  triangles  that  have  node  (i)  as  an  endpoint. 

Finally  the  new  locations  of  the  interior  nodes  are  calculated  by. 

xr'=xr+C’ 

Batina'^  states  that  the  predictor-corrector  procedure  is  more  efficient  than  performing  Jacobi 
iterations  alone.  This  is  due  to  a  fewer  number  of  iterations  required  to  achieve  acceptable 
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convergence.  Figure  3.8,  taken  from  Batina’s  paper,  shows  the  movement  of  the  grid  using  the 
dynamic  mesh  method. 


Initial  grid  Displaced  grid 

Figure  3.8  Illustration  of  the  dynamic  mesh  method  about  a  pitching  airfoil'^. 


This  concept  was  later  expanded  for  3-D  unstructured  meshes  by  Batina*''.  In  this  method  the 
mesh  is  modeled  by  tetrahedrons  with  springs  on  the  sides  as  shown  in  Figure  3.9. 


Figure  3.9  Tetrahedral  element  and  the  spring  equivalent. 


Robinson,  Batina,  and  Yang*’  modified  the  dynamic  mesh  to  work  with  structured  hexahedral 
grids  (Figure  3.10).  In  this  version  additional  springs  were  added  to  the  diagonal  of  each  cell  face  to 
control  cell  shearing.  A  factor  was  also  added  to  the  equation  for  the  spring  stiffness  for  the  ability  to 
control  the  spring  stiffness  of  the  cells  near  the  wing  to  prevent  distortion. 
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Singh,  Newman,  and  Baysaf  used  the  dynamic  mesh  algorithm  with  the  factor  for  spring  stifl&iess 
control  proposed  by  Robinson,  Batina,  and  Yang'®  on  a  2-D  unstructured  system.  In  this  paper  the 
number  of  elements  available  to  adapt  to  the  stmcture  were  limited  in  order  to  reduce  storage  and 
decrease  CPU  time.  The  authors  call  this  the  adaptive  window  procedure.  Figure  3.11  shows  an 
example  of  this  method. 


Initial  grid  Displaced  grid 

Figure  3.11  Unstructured  dynamic  mesh  about  a  pitching  airfoif 


Lesoinne  and  Farhat**  further  modified  the  dynamic  mesh  by  adding  fictitious  damping  and  mass 
to  the  spring  elements  to  create  a  moving  system  known  as  a  pseudo  structural  system.  The  authors  of 
this  paper  also  developed  a  method  for  analyzing  the  stability  of  the  dynamic  mesh  algorithm.  They 
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found  that  it  could  destabilize  the  fluid-structure  interaction  by  introducing  artificial  instabilities. 
Farhat”  added  to  the  pseudo-structural  system  by  placing  torsional  springs  at  the  vertex  of  the 
tetrahedron  to  prevent  a  vertex  from  interpenetrating  the  facet  of  a  tetrahedron  as  seen  in  Figure  3.12. 


An  advantage  of  the  dynamic  mesh  method  is  the  ability  to  specify  the  spring  stiffness  near  the 
structure.  This  gives  it  the  ability  to  handle  large  deflections  and  still  maintain  grid  quality.  One  of  the 
disadvantages  is  that  the  method  is  not  as  efficient  as  some  of  the  other  grid  updating  schemes.  This  is 
due  to  the  predictor  corrector  procedure  and  the  Jacobi  iterations.  The  method  also  becomes  more 
difficult  to  apply  and  less  efficient  as  extra  springs  are  added  to  the  network  to  control  cell  shearing  as 
discussed  in  Robinson,  Batina,  and  Yang”.  And  the  addition  of  the  torsional  springs  to  control 
interpenetration  as  discussed  by  Farhat”. 

3.2.3.4  Rigid  grid  method 

Tamura  et  al‘®  used  an  updating  procedure  implementing  ALE  by  making  the  structural  and  the 
fluid  grids  attached,  and  move  in  a  rigid  fashion.  This  method  is  similar  to  the  corotational  approachi* 
in  that  the  fluid  grid  is  attached  to  a  structure  undergoing  rigid  body  motion  and  the  entire  fluid- 
structure  system  moves  together.  What  makes  this  method  different  is  that  the  coordinate  axis  is  fixed 
in  space  and  does  not  move  with  the  computational  grid  as  seen  in  Figure  4-13. 

The  same  method  is  used  by  Selvam,  Govindaswamy,  and  Bosch”  to  model  the  incompressible 
flow  of  air  around  a  bridge  deck  section  to  determine  the  critical  velocities  of  the  structure.  Kandil  and 
Chung“  used  this  method  to  model  the  compressible  flow  of  air  aroxmd  an  airfoil. 
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In  this  method  the  entire  fluid  grid  is  rotated  about  an  origin  in  a  rigid  fashion  to  correspond  to  the 
changing  position  of  the  structure.  The  new  grid  is  found  by  simply  rotating  die  old  grid  in  a  rigid 
fashion  the  same  amount  of  angular  displacement  that  the  structure  was  moved.  Since  the  entire  fluid 


Figure  3.13  Movement  of  a  rigid  grid 


grid  is  moved  there  are  no  fixed  regions  (pure  Euler  coordinates),  therefore  the  whole  fluid  grid  is 
considered  ALE 

The  advantage  of  this  method  is  that  the  same  grid  is  used  throughout  the  entire  process.  Since  the 
grid  is  moved  in  a  rigid  fashion  the  element  areas  and  sizes  do  not  change  from  time  step  to  time  step. 
Therefore,  grid  quality  is  assiued  for  every  time  step.  The  method  is  also  easy  to  apply  due  to  the  grid 
displacement  being  the  same  as  the  structural  displacement.  No  numerical  methods  have  to  be 
introduced  to  predict  the  new  position  of  the  grid  making  the  method  computationally  efficient.  The 
disadvantage  that  the  outer  boundary  must  be  free  to  rotate  and  cannot  be  fixed  to  an  outer  boimdary. 


3.3  Theory  of  the  Deforming  Mesh 

A  gird  deformation  method'^’^’  which  is  suitable  for  aeroelastic  simulations  on  overset  grids  is 
used  in  the  following.  This  strategy  is  similar  to  TFI  in  that  it  is  algebraic  approach  based  on 
redefining  the  normal  grid  lines.  However,  unlike  TFI,  this  method  maintains  the  grid  quality  of  the 
initial  mesh  near  deforming  surfaces  imder  arbitrary,  moderate  deflections  and  rotations.  In  addition,  a 
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specified  region  in  the  far-field  may  be  held  fixed  so  that  the  grid  overlap  regions,  and  their 
connectivities,  remain  imchanged. 


3.3.1.  Calculation  of  the  Translations  and  Rotations 

Given  an  original  (starting)  grid  and  perturbed  sxuface  (assumed  to  be  at  k=l  here),  the  translation 
and  rotation  of  each  surface  node  can  be  computed  from  the  deflected  aerodynamic  surface. 

Defining  the  coordinates  of  the  original  and  perturbed  nodes  at  ? ,  ,k{=  )  as 


where  the  superscript  “1”  and/or  subscripts  “o”,  “p”  '•n'ote  t  =  l,  otigmal  and  present  respectively. 
The  translational  displacements  are  given  by 

AK'=X\-X[  (3-2) 

The  calculation  of  the  rotational  displacements  is  a  little  complex.  It  can  be  found  by  forming  an 
orthonormal  basis  for  die  original  surface  position  and  the  perturbed  surface  at  that  node.  The 
folloiving  four  steps  are  used  to  calculate  the  orthonormal  basis. 

(1)  Define  and  normalize  surface  vectors.  Suppose  the  given  node  is  1  shown  in  Figure  3.14.  The 

two  surface  vectors  and  b  are  defined  as 


Figure  3.14.  Definition  of  the  surface  vectors 


where  the  subscript  “2”,  “3”,  “4”,  and  “5”  denote  the  nodes  2, 3, 4,  and  5  respectively.  The  two  surface 


vectors  can  be  normalized  as 


0_  a 

d  —  II  11  5  D  II  11 


"H’  Wl 

If  the  present  surface  node  (i  )  is  a  boundary  node,  the  nodes  2,  3,  4,  and  5  defined  as  follows: 
(a)  set  incremental  numbers  of  nodes  2, 3, 4,  and  5  are  zeros,  i.e., 

4=0,  yi  =  0,(A:  =  2,3,4,5) 


(b)  Set 

i^=\,  ie 
i3  =  -l,  i^is 

js=-l,  s^js 

where  ie ,  i  ,  je ,  and  js  denote  the  ending  node  and  starting  node  in  the  i  and  directions 
respectively,  (c)  The  actual  nodal  number  of  1  through  5  are  given  by 

*  4  >  j  Jk ,  (^  =  42,3,4,5) 


Figixre  3.15.  Basis  on  the  given  surface 
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(2)  Form  and  unitize  the  surface  vectors.  As  shown  in  Figure  3.15,  the  summation  and  unitized 
vectors  of  the  two  surface  vectors  are  expressed  as 

+  (3-5) 

(3.6) 

“  hi 

(3)  Form  and  unitize  the  surface  normal  vector.  The  normal  vector  of  the  surface  is  the  cross- 
product  of  the  two  surface  vectors  a®  and  6“  which  are  shown  in  Figure. 

«,=a®x6®  (3.7) 


t'l  —  II  II 


(4)  Form  the  third  vector  for  the  basis 


The  orthonormal  basis  for  the  original  and  perturbed  surfaces  are  given  by 

^5  e,"i,  o-io 

Actually,  and  £  are  3x3  matrices.  The  corresponding  surface  rotation  matrix  can  be  defined  as 


3.3.2.  Update  a  Given  Grid  Line  (Biending) 

Each  normal  grid  line  is  moved  in  a  rigid-body  way  according  to  the  displacement  of  the  surface 
node  to  form  a  reference,  displaced  grid  line  defined  by 

X,  =  Xl+AX'+R  X^  -  Xl  (3.12) 

The  subscript  “r  ”  denotes  reference.  The  new  grid  line  is  constructed  by  blending  the  reference  grid 
line  and  the  old  grid  line.  The  blending  choice  is  arbitrary  but  is  best  done  in  arc-length  space  rather 
than  in  computational  space.  The  arc-length  for  each  node  is  defined 

=  (313) 

1=2 


where  5]  =  0 . 

A  cubic  blending  with  zero  slope  at  the  end  points  and  that  the  grid  transitions  smoothly  in  the  far- 


field.  This  can  be  written  as 
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V 


V  y 


J 


(3.14) 


where  represents  the  last  node  along  each  normal  line  that  is  allowed  to  deflect.  Fmally,  the  new 
position  of  each  grid  point  can  be  calculated  by  applying  the  blending  function  to  the  reference, 
displaced  grid  and  the  original  grid: 


(^  =  1,2,-) 


(3.15) 


3.4  Numerical  Simulations 

A  simply  supported  square  plate  is  considered.  The  side  lengths  and  thickness  are  non- 
dimensionalized  to  the  length  of  the  side.  The  detailed  information  of  the  plate  will  be  discussed  m 

pi 

Chapter  5.  The  plate  is  discretized  as  50x50  elements.  Possion  ratio  is  0.3.  ;t  =  -^-0.1. 

^  _  PjU_  _  gQ  jjjg  ijjtial  translational  velocities  on  all  the  nodes  are  0.01.  The  grid  coordinates 

D 

of  the  fluid  in  the  Z  direction  are  0.02  X/Z ,  where  IZ  is  the  integer  from  KS  to  KE.  The  other 
parameters  are  IE=JE=51,  IS=JS=KE=1,  NSF=1,  NFF=3,  and  IDEFM=3.  Two  cases,  KE=10  and  40 
are  considered.  The  moved  grids  for  the  two  cases  are  shown  in  Figures  3.16,  3.17,  and  3.18 
respectively. 

One  disadvantage  of  the  grid  deformation  approach  is  the  grid  will  overlap  when  the  deformation 
is  a  little  large  and  the  deformed  grid  in  the  boundary  is  not  usually  satisfied. 


3.5  Summary 

Different  methods  used  to  modify  the  fluid  equations  to  allow  the  movement  of  a  structure  within 
the  fluid  domain  are  described  in  this  chapter.  For  the  ALE  method,  various  methods  used  to  update  or 
move  the  fluid  grid  to  conform  to  a  changing  stmctural  boimdary  are  also  described  and  evaluated. 
After  that,  a  gird  deformation  method  which  is  suitable  for  aeroelastic  simulations  on  overset  grids  is 
presented.  This  strategy  is  similar  to  TFl  in  that  it  is  algebraic  approach  based  on  redefining  the 
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normal  grid  lines.  However,  unlike  TFI,  this  method  maintains  the  grid  quality  of  the  initial  mesh  near 
deforming  surfaces  under  arbitrary,  moderate  deflections  and  rotations.  In  addition,  a  specified  region 
in  the  far-field  may  be  held  fixed  so  that  the  grid  overlap  regions,  and  their  connectivities,  remain 
unchanged.  Numerical  simulation  shows  that  this  scheme  may  be  used  for  the  deforming  mesh  of  the 
panel  flutter. 
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(C) 

Figure  3.16.  The  original  grid  of  the  fluid  for  case  1:  (a)  in  3-D  coordinates; 
(b)  in  X-Y  coordinate  plane;  (c)  in  Y-Z  and  X-Z  coordinate  planes. 
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Figure 

(b)inX-Y 


(d) 

3.17.  The  deformed  grid  of  the  fluid  for  case  1:  (a)  in  3-D  coordinates; 
coordinate  plane;  (c)  in  Y-Z  coordinate  plane;  (d)  in  X-Z  coordinate  plane. 
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(b) 
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(c) 
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(d) 

Figure  3.18.  The  deformed  grid  of  the  fluid  for  case  2:  (a)  in  3-D  coordinates; 
(b)  in  Y-Z  coordinate  plane  for  1=25;  (c)  in  X-Z  coordinate  plane  for  J=25; 
(d)  in  X-Y  coordinate  plane  for  K=20. 
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CHAPTER  4 


COMPUTATIONS  OF  NAVIER-STOKES  EQUATIONS 
USING  FINITE  VOLUME  METHOD 


4.1  Introduction 

4.2  Governing  Equations 

4.3  Numerical  Procedures 

4.4  Boundary  Conditions  Implementation 

4.5  Computational  Results  and  Discussions 


4.1  Introduction 

The  fluid  dynamic  computation  must  have  the  fidelity  to  capture  the  relevant  flow  features  and 
provide  accurate  aerodynamic  loads  on  the  stmcture  in  developing  a  perfect  fluid-structure  interaction 
solver,  or  say  in  computational  aeroelasticity.  In  the  following  sections,  a  finite  volume  method  is  used 
to  establish  computational  solvers  for  Euler  equation  and  Navier-Stokes  equation. 


4.2  Governing  Equations 

4.2.1  Governing  Equation  in  Cartesian  Coordinates: 

4.2.1. 1  Non-Viscous  Flow 

The  governing  equations  of  non-viscous  fluid  flow  are  Euler  equations  and  the  corresponding 
continue,  energy  equations.  In  Cartesian  coordinate,  they  are  given  by 


U  F  G 
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>  ' 

'pu 

'pv 

pu 

,  F  = 

pu^  +  p 

puv 

pv 

puv 

pv  +p 

.PE. 

(E  +  p)u 

(E  +  p)v 

The  conservative  components  of  the  above  equations  are  given  by 

_£  +  JM+_1M  =  o 

dt  dx  dy 

ipu)  ^  (pm)  ^  (puv)  _  p 
dt  dx  dy  dx 

(pv)  ^  (puv)  ^  (pvv)  _  p 
dt  dx  dy  dy 

(pg)  ^  l(pE  +  p)u]  ^  (pE+  p)v]  _  Q 
dt  dx  dy 

where, 

g  =  e+-(M'+v') 
and  the  state  equation  is  defined  as 

p  =  (y-l)[pg--p(M"+v')] 


(4-2) 

(4-3) 

(44) 

(4-5) 

(4-6) 

(4-7) 

(4-8) 


4.2.1 .2  Viscous  Flow 

The  governing  equations  of  viscous  fluid  flow  are  Navier-Stokes  equation  and  the  corresponding 
continue,  energy  equations.  They  are  given  by 


U  F  G  . 
—  +  —  +  —  =  0 

dt  dx  dy 


(4-9) 


>  ’ 

pu 

pv  1 

pu 

,F  = 

pu^+p-r^ 

,  G  = 

1 

pv 

puv-T^ 

pv  +P-t:„ 

.pE. 

iE  +  p)u-ur„-v'C^+q,_ 

_(g  +  p)v-MT^-VT^+q^_ 

(4-10) 


The  conservative  components  of  the  above  equation  are  given  by 


_P+_(a)+_i^  =  0  (4-11) 

dt  dx 
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d(pu)  ^  djpm)  ^  d(puv)  ^  dp  ^  ^ 

^  dx  dx  dy 


(4-12) 

9(pv)  I  ^(P»v)  I  ^  I  I  (4-13) 

dt  etc  dy  dx  dy 

d(pE)  d[{pE  +  p)u]  d{pE+  p)v]  _  diuT^  ~^x)  ~^y)  ^4. J4) 

dt  ^  aje  dy  ~  dx  dy 

where,  according  to  the  Stokes  hypothesis, 

2 

and  A  is  second  viscosity,  for  gases  it  can  be  take  A  =  — p.  The  momentum  equations  are  given  by 
(ou\  (ouu^  (ouv)  p  2.M.  V  \T  .  r.,/  ^  fA  io\ 


(4-15) 

(4-16) 

(4-17) 


(pw)  , 

(puu)  1 

(PMV)  _ 

dt 

ar 

dy 

dx  dx 

(pv’)  , 

(P«v) 

(pvv)_ 

_^  +  _l 

at 

obc 

dy 

a^  a^ 

u  V. 


V  2  ,  M  .  V, 


For  convenience,  casting  the  governing  equation  (4-9)  in  nondimensional  form.  If  L  is  the 
characteristic  length,  and  other  characteristic  quantities  are  taken  to  be  ffeestream  values  denoted  by 
subscript  «>  ,we  may  define  the  dimensionless  variables,  denoted  by  an  asterisk,  as 


.  X  *  y  t 

"  ~T  '~l/ 


„♦=£,  V*=-,  /l*  =  ^  (4-20) 

L  L  L 

p-  =  p,  p-=-£-,  r=^,  e-  =  4 

^  L  p.v.  T.  Pi 

This  nondimensionalized  procedure  is  applied  to  the  compressible  Navier-Stokes  equation  (4-9),  the 
following  dimensionless  equation  are  obtained  as 


dt  dx 


(4-21) 
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Where  U\F\G' wt  the  vectors  given  by  equation  (4-10),  except  that  each  term  is  dimensionless, 
denoted  by  an  asterisk,  and  the  dimensionless  total  energy  per  unit  volume  is  given  by 


(4-22) 


The  components  of  the  shear  stress  tensor  and  the  heat  flux  vector  in  dimensionless  form  are  expressed 


*  2  3v* 

Re  3  etc*  3 

.  n’  du*  dv\ 

*  2  3m* 

^"'-ReSa/  3  3jc*'^ 


(7-1)M^  RePr  dx* 


=- 


(7-l)M^RePr  dy* 


(4-23) 


where  A/„,Re,Pr  are  the  freestream  Mach  niunber,  Reynolds  number  and  Prandlt  number 
respectively,  given  by 

=  ,  M.  Pr=^  (4-24) 


(4-24) 


jl  is  dimensionless  viscosity.  It  can  be  taken  as  constant  (=1.0)  if  the  temperature  variation  is  not 
very  large,  or  be  calculated  from  Sutherland’s  law 


.  ^(y.),.5  a+Q.L 

^  ^  ^  {T  +C,) 


(4-25) 


Where  r’is  dimensionless  temperature,  constant  C,  =  .  Here  T^is  freestream  temperature  (°K). 


The  perfect  gas  equation  of  state  become 
.  p*r* 

P  =^— r  >  ^  =■ 

yMl  1 


7(7 -1)M: 


(4-26) 


4.2.2  Governing  Equation  in  Non-orthogonal  Coordinates: 
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4.2.2.1  Non-Viscous  Flow 

From  the  Cartesian  equations,  we  can  get  the  governing  equations  in  non-orthogonal  coordinates: 

Continue  equation: 

P  I  tP") p  +  =0  (4-27 

3(  3^  3r|  3|  3>|  '' 

Because  -(Pi-Xj-,)!  =  [^^1  ■ «”«  '■‘k 

j  j  J  J 

3  p  3  (pl<)l.+(pvfe-,  3  ,(Pi'>l.+(pvX?,j_;,  (4,28 

3/7“^  J  dr]  J 

-(^)+^(^)+^(^)=o  (‘^29: 

dt  J  J  dr]  J 

Momentum  equations: 


then  we  have, 


3,^,  3(P!£^)^0 

dt  J  dq  J  dl]  J 


(4-28) 


(4-29) 


d{pu)  dipu^+p)^  d(pu^+P)„  ,  9(pwv).  Hpuv) 

IT*  3“^'  3n  ax  ' 


1  r^) + A  (ML!Ap\ + =  0 

dt^J^  d^^  J  dT]^  J 


(4-32) 


4.2.2.2  Viscous  Flow 

Momentum  equations: 

^  (y)  +  +  (|  (pIJ  .  +  +  (- 1  (pIJ,  +  ),  )v  +  ^,pVJ} 

+  ^{[puV  +  iUpr]J^  +  ipr]))u  +  i~(pi]J^+(jir]y)Jv+r]^p]/J} 

dr]  3  3 

=  ^(^[(j^/+^/)P«+|p^.<^.v]/J)  +  A(A[2(l4^77,+|^77,)At«+|p(^.t7,+l,^Jv]//) 
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- 1  )u  +  m.),  +  ^  +  ^yPV-^) 

dt  J  3  3 

+  ^{[pvF+((MT;,)„  -lipny)Ju  +  d(jiriy)y+(pri^)^)v+ii^p]/J} 
dri  3  3 


a  . 1 


+  ^(3-[(-nr  )PV  +  -W.J?,«)]/-0 
9tj  07]  3  '  3 


(4-34) 


Energy  equation: 

3  nP  s  8  ^  9  ,(«^xr +VT™ -?x),  9 

|-(f )+ |•(KP^ + .  ^([p^ =3j[ - — 1+  3^[ - j - J 

(4-35) 


4.3  Numerical  Procedures 


4.3.1  Finite  Volume  Method 

The  Governing  equations  (Navier-Stokes)  are  discretized  in  control  voliune: 


4.3.1 .1  Continue  equation: 

— (Pn‘-Pr;)  +  KP“)"l  -(P“)"'  ]A)'  +  [(Pv)”',-(pv)"MAi  =  0  (4-36) 

At  '+2->  '-i'i  ‘•■’*2  *’•'"2 

jc-direction  momentum  equation: 

At  ‘—’j  '•■'*2  '■^'2 


We  denote  the  node  (ij)  as  P  and  its  neighbor  nodes  (i+lj),  (i-lj),  (hj+l)  and  (iJ-V  W,N,S 
respectively,  and  the  surfaces  which  locate  at  (i-l-— ,7),  (/ — ,j),  (7»y  H  ).  ihj  )  between 


these  nodes  as  e,w,«,5  etc.,  as  shown  in  Figure  4.1. 
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Figure  4.1  Diagram  of  points  used  for  discretization 
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(4-42) 

The  ;c-direction  momentum  equation  is  rewritten  as 
ipur;'  -(p«)p 


At 


AxAy+[(puu),  -(puu)JAy  +  [(puv)„  -(puv),]Ax  =  -[p,iAy),  -p„(Av)«,]  + 

Re  3x  Re  dx  Re  dy  Re  dy 

IrP.  lAv  +  ir^r— )  -^(—)^A>c 

Upwind  scheme  is  used  for  the  convection  term,  and  central  scheme  is  used  for  diffusive  terms, 


then 


^^!^^AxAy-^^^^AxAy  +  xasx(F^Mpu)T'  -rmK(-F^,OXpu)f  -max(^;,0)(pM);*' 
At  At 


+maxi-FMpuyp'  +max(F„,0)(pM)r'  -max(-F„0)(p«);'^‘  -max(F„0)(pM) 
+  max(-F^,0)(pw)7'  =  -[Pei^y)e  “  p».(Ay)w]  + 

-[ipu)s-{pu)p]- .  -ipu)^]+ 


\«+l 


p^  Re(5r), 

pAy 

p„  Re(5r), 
Pe^y 

3p,  Re(&c), 
Pn^y 

p„  Re(&)„ 


•[(pw)w  -(P«)/>]- 


p„  Re(5r), 

Ps^y 
p,  ReCSr) 

pAy 


■[{pu)p-{pu)s]  + 


iipu), -ipu),]- ^ /^:^—[ipu)p -(pu)A  + 


3p„Re(5x)„ 


[ipv)„,-{pv)„J- 


Ps^y 


Ps  Re(5t), 


-lipv)se  -(Pv).„] 


(4^) 


Then,  the  x-direction  momentum  equation  in  discrete  form  is 

ap(pu)f  =fl£(pM)r'  +awip»yA  +«^,(p«);"' +«s(p«)r'  +*„  (4^5) 

where  the  coefficients  are  respectively 

Pe^y 


— 


Clw  — 


a,  = 


pAy 

P„  Re(5jc) 

PAy 


Pn 

PAy 

Ps  Re(&), 


+  inax(  -F^  ,0)  =  D^  +  max(  -F^  ,0) 
-+max(F,„0)  =Z)„  +max(F„,0) 

- + max(  -F„  ,0)  =  D„+  max(  -F„  ,0) 


+  max(  F^  ,0)  =D^+  niax(  F^  ,0) 
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inax(F,  ,0)  +  Z)„  +  inax(-F„  ,0)  +  +  max(F„  ,0)  +  £»,  +  max(-F,  ,0)  + 


AjcA>^ 

At 


Source  term 


K=(puyp^-[Pei^y)e-pA^y)J+b' 


Where  b  is  given  by 


b'  =  —E^^[{pu),-{pu),]-  -(P«)>r]+ 

3p,Re(5r),  3p„Re(5r)„ 


Pn  Re(&), 

The  iteration  equation  is 


p^  Re(&), 


=  (pM)7'  -a*  rest  at 


where  the  residual  is  given  by 


y-direction  momentum  equation: 

Similarly,  y-direction  momentum  equation  is  discretized  as: 

apipvy;'  =«£(pv)r'  +«»r(pv)r‘  +«v(pv)v'  +cis{pv)f  +b. 

Its  coefficients  are  all  the  same  as  that  in  x-direction  equation,  but  the  source  term  is 

AxAy 


fe.  =  (pv)^  ^  -  [P„  (Ax)„  -  p,  (Ax)  J  +  b^ 
A 


where, 


b'  = 


ttpv),  -(pv),] ^■'^,  ;-[(pv),-(pv),lt 


3p„  Re(^) 
Pe^y 


p,Re(^), 


[(P“)e«  -(pw)«]- 


3p,  Re(5y) 


p„Re(^),, 


[(P«)™  -(P“)».] 


(446) 

(447) 


(448) 


(449) 


res  =  aApu)f  -aAP^*)f  +«»r(P“)r'  +«v(P«r/  +«5(P«)r'  +*»  (4-50) 


(4-51) 


(4-52) 


(4-53) 


4.3.1 .2  Energy  equation: 

The  energy  equation  is 

^  a(p^ ^  a(pv^  ^  r aw  ^  3(£^] ^ 

3t  3x  dy  dx  dy  dx 

_ p _ rA(^)  +  i-(^)] 

(7-l)M^RePr  dx  dx  dy  ^ 


yx 


+  VT^)-H 
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(4-54) 


For  convenience  of  discretization,  rewrite  the  energy  equation  as 


^  ^  ^  3(£va  = -(^  +  +  i.  VT  )+ 1- (»I„  +  VI„  )  + 

dt  etc  dy  dx  dy  dx  oy 

MY  ^ - +  ^^^)]}  - 

RePr^ajc^9x  y(y-l)Afi  2  dy  dy  y(Y-l)Ml  2 

RePr  ax  ax  2  9^  9^  2 


(4-55) 


Because  specific  energy  E  =  ( - — - 5-+  “  ^  ^  )  >  ths  above  equation  is 

yfy-DM^  2 


7(y-l)M„  2 

^^a(pu£)^  a  a 

dt  dx  dy  dx  dy  dx  ay 

RePr‘ar^3^'  RePr*3^^9i'  2  Sy%'  2 


we  can  discrete  the  energy  equation  in  the  same  way  as  momentum  equation, 

OpipEy;'  =ap{pE)f  +a«,(p£:)^"'  +a^{pE)f  y  a^ipE)”;'  +h 
where  the  diffusion  conductance  in  the  coefficients  are  respectively 

n  D 

‘  p^RePr(5r)/  "  p^RePr(5r)„ 

n  i^«yAy  ^  ^sY^y 
"  p„RePr(&)/  *  p^RePr(5r)^ 


(4-56) 


(4-57) 


(4-58) 


Source  term  is 


h=b\yb^yybfyb^^yb\ 

bg  =  +  =  -[p,F,  -  +  p„F„  - p,F^  ] 

dx  dy 


(4-59) 


(4-60) 


(4-61) 
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+vt^)]AxA;^ 


bV  +vT,^)]AxAy 

dx 

Uu.Av  Adu  2  dv.  Hu„Ay  Adu  2  dv  flv.Ay  (du  dv  _  ^„Ay 

"iTV  Re  >  (4-63) 

8»  3v  4  2  8« 

Re  Re  dy  Re  S  3j^  3  3x  ”  Re  3  3}^  3  3jt  ■y 


^^u  ^  I  1  1  .  ^  1  /  M 

:v..Ax  A  1  2  1,  „  4  1  ,  ^21  , 


piv  Ax  A  \  2  1  1  .  2  1. 

”  (4-64) 

RePr  3x  dx  2  oy  oy  2 
_  jUy  Ay  ..^i+Vg  “p+Vp  ;uy  Ay  u],+v],  u^+Vy^ 
RePr(&:)^  2  2  RePr(&)^  2  2 

^y  Ax  ,ul+vl  «P+Vf-|  jUy  Ax  m^+v^  m.^+Vs-. 


- 

RePr  (Sy)„ 


^  -I  r'i  T'^P  ' 

2  RePr(fic)j  2 


bl^ipEf,^  <4-« 

In  the  above  discrete  equations,  the  terms  are  calculated  as  following: 

Ax  =  0.5(x£ ,  Ay  =  0.5(y;y -yj) , 

(A^)  +  ,  (Art. 


(Ax).  =0.x"'*  '"•'»).  (Ax),  =0.5(^^-^+^^^5-^) 

(&)^  —  X £  ^  Xp ,  (&)^  ^  Xp — x^r  >  (^)«  ”  >  (^X  ^  y p  y s 
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(5»:)„  =  (Ax)„  ,  (dx),  =  (Ax),  ,  =  {^y),  ,  =  (Ay), 

^ W  ~  Xsw 


(Ax),  +  .  (Ax),  =  0.5(- 


■) 


(Ay)„  =o.5(^^  y.K.+y^''  Ie^)  ,  (Ax),=o.5c^^  yw_^yEs 

0^  =O.5(0E+<t>p)  .0w=O-5(0p+0»f)  ><l>n  =0.5(0;^ +0;,)  ,  =O.5(0p+0s) 


(0)„.  =0.5( 

((l>)se  =0-5( 


,  ^NE 


-).W„w  =  0.5( 


0Ar 


) 


</>5  +  ^  ,  (^)  =  0.5(^^^!^  + 


(4-67) 


4.3.2  QUICK  Scheme 

The  accuracy  of  upwind  scheme  is  only  first-order  in  terms  of  Taylor  series  truncation  error.  The 
use  of  upwind  quantities  ensures  fliat  the  schemes  are  stable  and  obey  the  transportiveness  requirement 
but  the  first-order  accuracy  makes  them  prone  to  numerical  diffusion  errors.  Such  errors  can  be 
minimized  by  employing  higher  order  discretisation. 

The  quadratic  upwind  interpolation  for  convective  kinetics  (QUICK)  scheme  uses  a  three-point 
upstream- weighted  quadratic  interpolation  for  cell  face  values.  The  face  value  of  0  is  obtained  from  a 

quadratic  function  passing  through  two  bracketing  nodes  (on  each  side  of  the  face)  and  a  node  on  the 
upstream  side.  The  Hayase  QUICK  scheme  can  be  summarized  as  follows: 
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+-(30?  “ 2^j  - (pss) 

for  F,  >  0 

=<l>P+-0<l>N  -'^P  -^s) 

for  F„  >  0 

~^P  ~'^^P  ~*Pn) 

for  <  0 

~^N  +“(3^/>  ~^Nn) 

for  F„  <  0 

(4-68) 

The  QUICK  scheme  has  greater  formal  accuracy  and  retains  the  upwind  weighted  characteristics. 
The  resultant  false  diffusion  is  small  and  solutions  achieved  with  coarse  grids  are  considerably  more 
accurate  than  those  of  the  upwind  or  hybrid  schemes.  The  QUICK  scheme  can,  however,  give  (minor) 
undershoots  and  overshoots  as  is  evident  in  computational  results  (Figure  4.5). 

4.4  Boundary  Conditions  implementation 


4.4.1  Inlet 

For  the  non-dimensional  equations,  at  inlet  the  following  conditions  are  given  by 
M  =  1.0  ;  v=0.0  ;  p  =  •,P  =  P^,T=T„ 

^  1  ■  j._  1  1 

ymJ  ’  yir-Wl  2 


(4-69) 


4.4.2  Onlet 

In  flow  direction,  gradients  of  variables  are  zero,  which  is; 

-1  =  0  (4-70) 

I 

4.4.3  Solid  Wall 

4.4.3.1  Euler  equation  boundary  conditions 

The  basis  is  impermeable  boundary  condition: 

F  =0  (4-71) 

n 

If  the  coordinates  are  non-orthogonal,  we  have. 
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F„=F-n 

=  ««  +  (^4  ■  S )  +  •  ^« )  +  ) 

+  Fj/J^(ef  •e|)  +  Kjn,(ef  -e^) 

=  F^n^+F^n^  +  F^n^ 

+  (F^n^  +F'^«^)[cos(|,7?)]  +  (J^««f  +ff«4)[cos(^,0]  +  (^n"f  +^f«t,)[cos(r?,0] 

=  +F5«f  [cos(^,C)]+f;”c[cos(t7,C)] 

(4-72) 

Here,/!  is  the  surface  normal  direction,  n  =  n^e^  +njje^  +  ^  rf  Coordinates  are  on  the  solid 

surface. 

For  two-dimensional  flow  problem,  the  above  equation  can  be  write  as: 

F„  =  F^n^  +  F^  n,  [cos(^ ,/))]  (4-73) 

Here  the  rj  coordinate  is  on  the  solid  surface. 


Planar  waU: 

The  I  t]  components  of  velocity  .F^jand  other  parameters  are  extended  to  boimdary  by  linear 
interpolation  or  symmetry  techniques  are  used.  F^  is  obtained  with  the  above  equation.  For  example, 
when  synunetry  techniques  are  used,  the  boundary  conditions  are  as  follows: 


P-i  =  F+i  .  P-i  =  P+i 


(KU  =(n).. 


iV(U=iVi\, 


(4-74) 

(4-75) 


(Kh  =-(K).,  ^ 

{Fj/if  +F^/ij[cos(^,^)]  +  Fq/if  [cos(7),^)]}_,  =-{F^n^  +F^n^[cos{^,C)]  +  F^n^[cos{riX)]}+\ 

(4-76) 

Here  point  -1  is  the  imaginary  point  outside  the  boundary. 


Curvature  wall: 

For  Curvature  wall,  the  pressure  boundary  condition  should  be  corrected  with  the  following 
equation: 


(4-77) 


Here  R  is  radius  of  curvature.  Then  the  pressure  condition  is  corrected  as: 


110 


CHAFER  4:  Computations  for  Navier-Stokes  Equations  Using  Finite  Voltime  Method 


Here  A«,  represent  the  distance  between  points  1  and  -1.  The  density  boundary  condition  is: 


(4-78) 


(4-79) 


4.4.3.2  Navier-Stokes  Equation  Boundary  Conditions 

On  the  solid  wall  boundary,  the  viscous  flow  no-slip  conditions  are  applied. 

u=0  ,  v=0  (4-80) 

The  temperature  boimdary  condition  is  adiabatic  wall  condition  or  given  the  freestream  temperature 

7L 


—  =  0  or  T  =  T^ 


(4-81) 


4.5  Computational  Results  and  Discussions 

The  two  numerical  methods,  Euler  solver  and  Navier-Stokes  solver,  described  in  the  previous 
section  were  combined  into  a  single  code  in  such  a  way  that  we  could  readily  select  between  one 
method  and  the  other  in  solving  a  specified  flow  problem. 

The  problem  of  the  flow  in  a  channel  with  a  circular  arc  “biunp”  was  chosen  to  evaluate  the  code 
for  subsonic,  transonic,  and  supersonic  steady  state  modeling.  This  particular  problem  is  well  suited 
for  code  development  and  testing.  The  geometry  and  the  grids  are  easy  to  generate  accurately  and  the 
problem  symmetry  and  geometrical  simplicity  aid  the  interpretation  of  the  results.  Two  circular  arc 
bump  thickness-to-chord  ratios  were  used:  10%  for  subsonic  and  transonic  modeling,  4%  for  the 
supersonic  model.  To  the  upper  and  lower  boundaries  of  the  channel,  the  solid  wall  boundary 
conditions  were  applied.  The  inflow  boundary  is  on  the  left  side. 

Figure  1  is  the  computational  geometry  for  the  flow  model.  Figure  2  shows  isomach  lines  of  the 

steady  flow  solution  for  an  upstream  Mach  number,  =  0.5  ,  obtained  using  the  Euler  solver.  The 
isomach  lines  distributions  are  nearly  symmetry.  That  is  the  verification  of  correctness  of  the 
computational  code.  The  contour  is  not  smooth  enough,  that  may  be  result  from  the  hyperbolic 
characteristics  of  Euler  equation  and  the  first-order  boundary  condition  implementation.  Improving  the 
boimdary  conditions  to  second-order  or  adopting  finer  grids  can  make  it  better.  Figure  3  shows  the 
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channel  upper  and  lower  surfaces  Mach  number  distribution  for  the  same  test  case.  A  little  asymmetry 
is  due  to  the  influence  of  the  sharp  leading  and  trailing  edge.  Figure  4  and  5  present  the  contour  of 
Mach  number  and  the  surface  Mach  number  distribution  of  the  transonic  flow  solution  for  the  same 
10%  thick  arc  bump  for  M_„  =0.675  using  the  Euler  solver. 

The  computational  code  can  simulate  the  sharp  shock  profile  very  well  that  appears  on  the  bump, 
except  that  there  is  a  small  spurious  oscillating.  That  is  the  dispersion  introduced  by  QUICK  scheme 
along  with  its  improvement  in  accuracy.  Figure  6  and  7  present  results  for  the  supersonic  flow  in  a 
channel  with  a  4%  thick  arc  bump  for  =  1.65  obtained  with  Euler  code.  The  shock  wave,  their 
reflection  and  interaction  can  be  simulated  pretty  well.  The  quality  is  improved  when  the  grids  get 
finer.  Figure  8  and  9  show  the  contour  of  Mach  number  and  near  surface  Mach  number  distribution  of 

subsonic  flow  solution  for  the  10%  thick  arc  bump  for  upstream  Mach  number  M_„  =  0.5  using 
Navier-Stokes  equations. 

The  Mach  number  contour  is  asymmetry,  not  as  that  of  results  of  Euler  solver.  The  reason  is  that 
viscous  boundary  layer  developing  along  the  surface  and  that  the  boundary  separation  near  the  trailing 
edge  of  the  bump  (Fig.8).  The  boundary  layer  and  its  separation  also  change  the  flow  passage  in  non- 
viscous  sense,  therefore,  the  Mach  munber  increases  and  the  flow  pattern  is  not  symmetry  anymore. 
Figure  10  and  11  present  the  results  for  transonic  flow  in  the  chaimel  for  the  same  test  case  with 

M_^  -  0.675  using  Navier-Stokes  equation. 

The  shock  wave  is  a  little  smeared  because  of  the  physical  viscosity,  and  at  the  same  time  the 
spurious  oscillating  does  not  appear  also  due  to  the  physical  viscosity.  There  is  a  shock  and  boundary 
interaction  downstream  of  the  bump.  The  boundary  layer  and  its  interaction  with  the  shock  affect  the 
flow  field  in  the  channel.  Figure  12-15  shows  supersonic  flow  over  a  flat  plate  with  0%  thickness  for 
upstream  Mach  M =  1.8  using  Navier-Stokes  equation.  The  upper  boundary  is  free  stream. 

A  shock  wave  spears  in  the  middle  of  the  plate,  its  action  with  the  boundary  layer  generates  an 
expansion  shock  downstream  of  it.  The  two  shocks  dominate  the  flow  field  pattern  over  the  flat  plate. 
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Figure  4.2  The  computational  geometry  of  flow  over  a  bump 
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Figure  4.3  Contour  of  Mach  number  of  flow  over  a  bump 
(Inflow  Mach  number  M_^  =  0.5 ,  Euler  solution) 
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Figure  4.4  Mach  number  distribution  on  the  upper  and  lower  surface 
(Inflow  Mach  number  M_„  =  0.5 ,  Euler  solution) 


Figure  4.5  Contour  of  Mach  number  of  flow  over  a  bump 
(Inflow  Mach  number  =  0.675 ,  Euler  solution) 
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Figure  4.6  Mach  number  distribution  on  the  upper  and  lower  surface 
(Inflow  Mach  number  M_„  =  0.675 ,  Euler  solution) 


X 

Figure  4.7  Contour  of  Mach  number  of  flow  over  a  bump 
(Inflow  Mach  number  =  1.65 ,  Euler  solution) 
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Figure  4.8  Mach  number 
(Inflow  Mach  nw 
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Figure  4.9  Contour  of  Mach  number  of  flow  over  a  bximp 
(Inflow  Mach  number  M_^  =  0.5 ,  Navier-Stokes  solution) 
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Figure  4.10  Mach  number  distribution  near  the  upper  and  lower  surface 
(Inflow  Mach  number  =  0.5 ,  Navier-Stokes  solution) 


Figure  4. 1 1  Contour  of  Mach  number  of  flow  over  a  bump 
(Inflow  Mach  number  M_„  =  0.675 ,  Navier-Stokes  solution) 
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Figure  4.12  Mach  number  distribution  near  the  upper  and  lower  surface 
(Inflow  Mach  number  M_„  =  0.675 ,  Navier-Stokes  solution) 
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1.5  1 
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N  over  flat  plate  (Inflow  mach  number  M_^  =3.0) 


V  over  flat  plate  (Inflow  mach  number  M_^  =3.0) 
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CHAPTER  5 


THREE  DIMENSIONAL  AEROELASTIC  SOLVER 
FOR  NONLINEAR  PANEL  FLUTTER 


5.1  Nomenclature 

5.2  Introduction 

5.3  Aerodynamic  Theories 

5.4  Governing  Equations  of  Thin  Plate 

5.5  Numerical  Procedure 

5.6  Results  and  Discussions 

5.7  Summary 


5.1  Nomenclature 

R  =  membrane  rigidity  of  plate  defined  in  equation  (5-33); 

c  =  speed  of  sound  in  air; 

D  =  bending  rigidity  of  plate  defined  in  equation  (5-33); 

E  -  total  specific  energy  per  unit  volume  defined  in  equation  (5-22); 

=  elastic  modulus  of  plate; 
h  -  thickness  of  plate; 

H(^l  ~  source  vector  term; 

J  =  Jacobian  of  the  transformation; 

/  =  side  length  of  plate;  characteristic  length; 

=  free-stream  Mach  number; 

=  dimensionless  air  pressure; 

=  Prandtl  number,  0.73  for  standard  air; 
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=  aerodynamic  pressure,  positive  in  the  direction  opposite  to  w ; 

=  aerodynamic  pressme  due  to  plate  motion; 

=  ‘fextemal”  aerod5mamic  pressure  independent  of  plate  motion,  e.g.  turbulent  pressure 
fluctuations; 

=  Reynolds  number,  p^VJ  j ; 

=  dimensionless  time; 

=  nonlinear  internal  force  coefficient  defined  in  equation  (5-38); 

=  dimensionless  static  temperature; 

=  dimensionless  deformations  of  the  plate  in  x  z  directions; 

=  dimensional  deformations  of  the  plate  in  x  y  z  directions; 

=  dimensionless  displacement  vector  of  the  plate  in  x  direction:  u-  m,  ,1/2,"  •  ,«„  ; 

=  velocity,  V  =  vj  +  v^j  +  vjc ; 

=  free-stream  velocity; 

=  three  components  of  the  velocity; 

=  dimensionless  displacement  vector  of  the  plate  in  direction:  v=  v,,V2,"',v„  ; 

=  dimensionless  displacement  vector  of  the  plate  in  z  direction:  w  = 

=  dimensionless  Cartesian  coordinates; 

=  constants  for  Newmark-  p  integration  method; 

=  Kronecker  delta  function; 

=  velocity  potential  function; 

=  ratio  of  specific  heat,  1.4  for  perfect  gas; 

=  dimensionless  dynamic  pressure  defined  in  equation  (5-38); 

==  molecular  viscosity  coefficient; 

=  mass  ratio  defined  in  equation  (5-38); 

=  dimensionless  mass  density  of  air; 

=  dimensionless  mass  density  of  plate; 

=  Poisson’s  ratio; 
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T  =  time  in  computational  space; 

^  JJ  ^  =  Coordinates  in  the  computational  space; 

Subscripts: 

=  air; 

b  =  boundary; 

s  =  stmcture  (plate); 

=  viscous; 

oo  =  free-stream  value; 

Overhead: 

^  =  denote  the  dimensional  parameter 


5.2  Introduction 

An  important  goal  of  computational  aeroelasticity  is  to  impact  the  design  process  with  simulations 
of  full  aircraft  configurations.  In  developing  a  fluid-structure  interaction  solver,  three  sets  of  governing 
equations,  as  well  as  their  coupling,  must  be  considered.  These  are  the  fluid  dynamic,  structural 
dynamic,  and  fluid-structure  interface  equations.  The  fluid  dynamic  equations  must  have  the  fidelity  to 
capture  the  relevant  flow  features  and  provide  accurate  loads  on  the  structure.  The  structural  dynamic 
equations  must  in  turn  be  capable  of  modeling  the  structural  deformation  properties  under  a  given 
time-varying  load.  Since,  in  general,  the  aerodynamic  and  the  structural  discretizations  are  not 
identical,  accurate  and  stable  structural  interface  equations  must  be  established  to  transfer  the 
aerodynamic  loads  to  the  structure  and  the  structural  deformations  to  the  aerodynamic  mesh. 

When  considered  independently,  high-fidelity  fluid  dynamic  and  structural  dynamic  solution 
algorithms  have  become  mature  for  flows  or  structural  members  exhibiting  nonlinearities.  However, 
application  of  these  solvers  in  a  time  accurate,  multidisciplinary  environment  requires  further  study  in 
order  to  elucidate  various  issues  arising  in  the  coupling  of  these  high-fidelity  solvers. 

Historically,  researchers  interested  in  dynamic  aeroelastic  computations  have  taken  well-validated, 
implicit  Navier-Stokes  algorithms  developed  to  solve  complex  flows  over  three  dimensional,  rigid 
bodies,  and  extended  them  to  include  aeroelastic  effects.  The  most  common  practice  is  to  simply  lag 
the  effects  of  moving/deforming  structures  by  one  time  step^’^  allowing  current  algorithms  to  be  used 
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in  updating  the  aerodynamic  variables.  Bendiksen  and  Hwang^  pointed  out  that  when  took  this 
approach  unknown  phrase  and  integration  errors  were  introduced  leading  in  some  cases  to  incorrect 
prediction  of  the  stability  behavior  of  the  fluid  structure  interaction  system.  To  overcome  this  problem 
Bendiksen  and  Davis'*  took  the  alternate  approach  to  develop  a  new  tightly  integrated  algorithm  in 
which  the  fluid  and  structures  are  modeled  as  a  single  dynamical  system.  Although  this  approach  has 
been  shown  to  eliminate  the  lagging  errors,  it  requires  the  development  of  an  entirely  new  solver.  This 
would  entail  large  additional  investments  by  the  aircraft  industry  for  developing  and  validating  new 
flow  solvers,  and  would  result  as  well  in  the  loss  of  the  significant  user  experience  already 
accumulated  with  existing  aerodynamic  and  structural  tools. 

An  attractive  alternative  method  for  eliminating  these  phrase  and  integration  errors,  while  utilizing 
existing  fluid  dynamic  and  structure  dynamic  algorithms,  is  implementing  Newton-like 
subiterations^’^.  Subiterations  can  eliminate  the  errors  from  the  linearization  and  factorization,  as  well 
as  from  the  lagging  boundary  conditions  and  turbulence  models.  The  added  computational  cost  of 
subiterations  is  typically  an  additional  solution  vector  storage,  and  each  subiteration  is  equivalent  in 
workload  to  a  time  step  of  the  baseline  algorithm.  The  result  is  a  fiilly  implicit  coupling  between  the 
fluid  and  structures  without  having  to  develop  a  completely  new  tightly  coupled  solver. 

Morton,  Melville,  Visbaf  applied  this  approach  to  an  elastically  mounted  cylinder  in  a  uniform 
crossflow.  The  subiteration  methodology  was  shown  to  reduce  the  structural  coupling  errors  and 
allowed  higher  order  accurate  time  integration  scheme  to  be  used  with  relatively  minor  changes  to  the 
baseline  aerodynamic  solver.  The  overall  workload  of  the  algorithm  was  improved  by  a  factor  of  ten 
over  traditional  first  order  lagged  approaches  due  to  demonstrated  second  order,  coupled,  temporal 
integration.  The  structural  dynamics  solver  was  comprised  two  linear,  second  order  differential 
equations  describing  motion  in  the  horizontal  and  vertical  axes  and  were  integrated  in  time  with  the 
same  method  used  in  the  aerodynamic  code.  The  fluid-structural  interface  equations  were  trivial  for 
this  simple  model.  Melville,  Morton,  and  Rizzetta^  also  have  used  this  technique  to  couple  a  three 
dimensional  Navier-Stokes  code  with  a  general,  linear  second-order  structural  model.  This  solver  has 
been  applied  successfully  to  the  problems  of  transonic  wing  flutter^,  tail-buffet^,  and  limit-cycle 
oscillations*^. 

In  order  to  identify  pertinent  issues  related  to  coupling  nonlinear  structural  dynamic  and 
aerodynamic  equations,  a  model  problem  that  is  geometrically  simple,  has  reasonable  computational 
requirements,  yet  requires  high-fidelity  fluid  and  structural  solvers  is  desired.  One  such  model 
problem  is  nonlinear  panel  flutter. 
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Only  a  few  computational  studies  have  recently  been  considered  the  nonlinear  flow  effetcs  in 
panel  flutter.  Davis  and  Bendiksen",  Davis Bendiksen  and  Davis'*  have  employed  an  improved 
modeling  of  the  aerodynamics  by  tightly  coupling  the  Euler  equations  with  a  nonlinear  finite  element 
model  for  two-dimensional,  Transonic  panel  flutter  problems.  Selvam,  Visbal,  and  Morton‘S  extended 
the  subiteration  methodology  to  more  complex  structural  models  including  nonlinear  effects  and 
requiring  the  coupling  of  different  time-integration  schemes  for  the  fluid  and  the  structure.  In  addition, 
they  provided  some  insight  into  the  effects  of  viscosity  on  thin  panel  flutter  in  the  transonic  regime. 
Gordnier  and  Visbaf'*  extended  the  preliminary  work  of  Selvam,  Visbal  and  Morton  for  two- 
dimensional  panel  flutter.  Viscous  effect  on  the  flutter  dynamic  pressure  for  supersonic  and  subsonic 
flow  was  discussed.  Additionally,  some  computations  were  performed  for  inviscid,  three-dimensional 
panel  flutter. 

Several  commonly  used  approximate  aerodynamic  theories,  linearized  potential  flow  theory,  qusi- 
steady  supersonic  theory,  first  order  piston  theory  and  full  three  order  piston  theory,  are  first  listed  in 
Section  3  for  comparison  purposes.  The  valid  range  of  Mach  number  is  provided  for  each  theory.  The 
unsteady,  compressible,  three-dimensional  Euler  equations  are  also  provided  in  this  section  to  predict 
the  aerodynamic  pressure  acting  on  the  panel  due  to  its  deformation.  To  reduce  the  interaction  error 
between  the  fluid  and  structure,  the  expression  for  the  Geometric  Conservation  Law  (GCL)  is 
discussed.  In  Section  4,  the  governing  equations  based  on  the  von  Karman  large  deflection  theory  are 
provided.  All  these  governing  equations  are  nondimensionalized  by  the  appropriate  combination  of 
freestream  density,  velocity  and  the  length  of  panel.  In  Section  5,  the  Beam- Warming,  alternate- 
direction,  implicit  scheme  is  listed  to  solve  the  flow  equations.  The  finite  difference  method  and  the 
Newmark-  integration  scheme  are  used  to  solve  the  governing  equations  of  nonlinear  panel.  The 

Newton-like  subiteration  is  also  implemented  in  this  section  to  eliminate  the  lagging  errors  associated 
with  the  exchange  of  the  pressure  and  deformations  between  the  fluid  and  panel  at  their  interface.  The 
nonlinear  panel  flutter  under  Euler  flow  is  simulated  for  a  square  panel  in  Section  6.  The  amplitudes  of 
the  l  imit  Cycle  Oscillation  (LCO)  are  compared  with  those  obtained  from  full  linearized  potential 
flow  theory  and  qusi-steady  linear  piston  theory.  The  flutter  frequencies  are  also  discussed  and 
compared  with  their  natural  frequencies. 


5.3  Aerodynamic  Theories 

5.3.1  Approximate  Aerodynamic  Theories 
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Normally  the  aerodynamic  pressure  may  be  considered  as  the  sum  of  two  parts,  one  given  by  the 
pressure  fluctuations  on  the  plate  in  the  absence  of  any  plate  motion,  e.g.,  due  to  turbulent  boimdary 
layer  fluctuations,  and  the  other  due  to  the  plate  motion  itself”.  The  superposition  of  these  two 
contributions  to  form  the  total  aerodynamic  pressure  implies  that  the  plate  motion  and  the  consequent 
portion  of  the  aerodynamic  pressure  are  sufficiently  small  so  that  the  turbulent  pressure  fluctuations 
themselves  are  not  substantially  modified.  Hence, 

#  =  4p"+4p^  (5-1) 

Ap^  is  taken  as  prescribed  from  experiment  or  separate  analysis  while  Ap^  may  be  related  to  the 
deformation  of  the  plate.  The  fiill  expression  of  the  Ap^  was  defined  by  Dowell”. 

As  disclosed  by  the  survey  papers,  a  vast  quantity  of  literature  exists  on  panel  flutter  using 
approximate  aerodynamic  theories  without  considering  the  effect  of  Ap^ .  Several  frequently  used 
theories  are  listed  in  the  following.  Their  valid  ranges  of  Mach  number  were  summarized  by  Mef®. 


5.3.1.1  Linearized  Potential  Flow  Theory 

For  air  flow  with  Mach  number  close  to  one  and  less  than  five  (1  <  Af„  <  5 ),  the  full  linearized 
inviscid  potential  theory  aerodynamics  is  usually  employed”.  The  aerodynamic  pressure  loading  is 
given  by 


fd(l> 


+y. 


d(l> 

dx 


where  the  velocity  potential  flmction  0  must  satisfy 

subject  to  the  boundary  conditions 


dip 

dz 


li=0 


dw  ^  dw 
0 


on  plate 
off  plate 


(5-2) 


(5-3) 


(5-4) 


5.3.1. 2  Qusi-Steady  Supersonic  Theory 

The  aerodynamic  theory  employed  for  most  part  of  the  panel  flutter  at  higher  supersonic  Mach 
munber  (A/^>1.6)  is  the  quasi-steady  first  order  piston  aerodynamic  theory  by  Ashley  and 
Zartarian”.  The  aerod5mamic  pressure  as  given  in  this  theory  is  given  by 
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2q^  ^dw  MI-2  1  dw^ 

P  3Jc  Ml  -1  V  dt 
^  V  -  «o  y 


(5-5) 


where  qa~  stream  dynamic  pressure;  If  the  aerodynamic 

damping  is  neglected  in  Equation  (5-5),  the  quasi-static  Ackeret  aerodynamic  theory,  also  know  as  the 
static  strip  theory,  is  simplified  to 


2qa  w 
■  p  =  — — — 

p  dx 


(5-6) 


5.3.1. 3  First  Order  Piston  Aerodynamic  Theory 

The  simplified  formula  of  the  aerodynamic  pressure  within  the  first  order  piston  theory  has  the 
expression'^’*’ 


p-p 

M_ 


^dw  1 

dx  V  di 


(5-7) 


This  approximation  is  usually  used  for  the  supersonic  and  hypersonic  flow  (2.5  <  M„  <5  by  Krause 


and  Dinkier^*,  or  V2  <  <  5  by  Mei'*). 


5.3.1. 4  Full  Third  Order  Piston  Aerodynamic  Theory 

In  the  hypersonic  regime  ( >5),  the  unsteady  full  third-order  piston  theory  aerodynamics'®  is 


used  to  develop  the  aerodynamic  pressme  given  by 


^  2l 


1  dw 
dx  dt 


^dw  1  dw^ 

_ L _ 

,  iy+ml 

^dw  1 

3" 

4 

^dx  V^di  ^ 

*  12 

dx  V  dt 

K  °°  J 

where  v]  is  the  ratio  of  specific  heat. 


(5-8) 


5.3.2  Aerodynamic  Theory  Based  on  Navier-Stokes  Equations 

The  aerodynamic  governing  equations  are  the  unsteady,  compressible,  three-dimensional  Navier- 
Stokes  equations  written  in  nondimensional,  strong-conservation  law  form  employing  a  general  time 
dependent  transformation 


127 


CHAPTER  5:  Three  Dimensional  Aeroelastic  Solver  for  Nonlinear  Panel  Flutter 


n-r^x,y,z,i) 

C  =!^{x,y,z,t) 
t  =  t 

The  resulting  system  of  governing  equations  is  repressed  as 

du  d(-  1  -"I  ^  hVh  r5io^ 

dt  Re,  J  dr][  Re,  j  J 

The  source  vector  term,  Hqcl  >  defined  as 


Hgcl^^U  +f^l  +f^ 

®  dx  J  L  /  I  J 


(5-11) 


is  applied  in  order  to  account  for  errors  associated  with  the  GCL  and  will  be  discussed  later.  J  is 
Jacobian  of  the  transformation^*. 

Vector  quantities  appearing  in  Equation  (5-10)  are  defined  as 


U=-U 

J 


(5-12) 


^  =  y  (n,/  +  J7.X^  +  +  n,z^) 


(5-13) 


(5-14) 


#  =  i(C,tJ  +  CF  +  C.,<?  +  C,z^) 

G,=^^^Fy+n.yG,-^%H,) 

With  this  formulation,  the  vector  of  dependent  variables  U  is  given  as 

Po''^  Pa^y  P-'^^  P-^"^ 


(5-15) 


(5-18) 


(5-19) 


and  the  flux  vectors  are 
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PaK  P  Pa^x^y  Pa^x^z 

F  =  Pa^x^y  »  ^  =  Pa^l  P  >  H  =  Pa^y^z 
PaVx^z  Pa^y^z  Pa^l  +  P 

(P^E  +  px\  [iPaE  +  P>y\  [(p,E+p)v^ 

=  0  Txy  T„  +V^T^  +V,T„  -qx  ’’ 

=  0  T^xy  -^yy  '^yz  xy  + 

^v=0  T„  V,T„+VT  +V,T„-^/ 


(5-20) 


(5-21a) 
(5-2  lb) 


(5-2 Ic) 


where  the  total  specific  energy  per  unit  is  defined  as 
T  (v^  +  +  v^) 

E  = - - - r  +  ^  (5-22) 

y(y  -  1)M  £  2 

All  variables  have  been  normalized  by  the  appropriate  combination  of  freestream  density,  velocity  and 
a  characteristic  length,  that  is. 


X  y 

X-.,  y-..  z-.,  . 


(5-23a) 


Vx  ^y  p 

Vx=-^  ,  V„  =  -/-  ,  V=-!^,  U  = 

K.  ^  K. 


(5-23b) 


-  _A  --  P  T-  — 

p-  pfy  t 

Components  of  the  stress  tensor  and  the  heat  flux  vector  may  be  expressed  as 
^  du,  5m  /  2  j.  dux  ^ 


_J _ Tp 


(5-23c) 


(5-24) 


where  i 


;  Mi,M2,M3  =  v^,v  ,v^  and  x,,X2,*3  =x,y,z.  Sutherland’s  law  for  the  molecular 


viscosity  coefficient  p  and  the  perfect  gas  relationship 

pT 

P=-^ 


are  also  employed,  and  Stokes’  hypothesis  for  the  bulk  viscosity  coefficient  is  assumed. 
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5.3.3  Geometric  Conservation  Law 

The  relationship  between  the  governing  equations  and  the  GCL  is  discussed  in  the  following.  The 
nondimensional  Cartesian  governing  equations  can  be  expressed  as 


— +  —  F — +  — 
3x1  Rzi  ”  I  3^ 


1  \  3  f  1  1 

+f  H-  —  H,  =0 
Re,  I  3zl  Re, 


(5-27) 


Using  chain-rule  differentiation, 

_o  =  _o + _Q^+ _QJI + _o  _£ 

3r  3t  3|  dt  djj  dt  3C  3t 

_Q^_0^+_0jL+_0_i 

3x  3^  etc  dr]  dx  3^  3x 

_Q  =  _0_i+_0JL+Ji_£ 

dy  d4  ^  dr]  dy  d^ 


(5-28a) 

(5-28b) 

(5-28c) 


_Q=_Q_i+_OjL+_Q_£ 

dz  d^  dz  dr]  dz  d^  dz 


(5-28d) 


and  premultiplying  by  the  inverse  of  the  transformation  Jacobian,  J ,  equation  (5-27)  becomes 


3t  ^3| 


All  four  terms  on  the  right  hand  side  of  equation  (5-29)  vanish  analytically.  The  difficulty  arises 
when  discrete  representations  of  the  temporal  and  spatial  derivatives  are  used.  The  discrete  form  of  the 
last  three  terms  are  zero  when  central  differences  are  used  for  all  metric  calculations  in  3-D. 
Unfortunately,  this  is  not  true  for  the  first  expression  due  to  the  mixed  temporal  and  spatial 
derivatives.  The  first  term  set  to  zero  is  referred  to  in  the  literature  as  the  The  most 
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straightforward  approach  of  accounting  for  the  GCL  is  to  simply  include  this  term  in  the  discrete 
governing  equations. 


5.4  Governing  Equations  of  Thin  Plate 


The  governing  equations  of  the  thin  plate  are  the  well  known  von  Karman  equations  for  large 
deflection.  Derivations  of  these  equations  may  be  found  in  a  number  of  sources^'*.  For  the  von  Karman 
theory  the  plate  is  assumed  to  be  isotropic,  of  uniform  small  thickness  and  initially  flat.  The  normal 
deflection  of  the  plate  is  assumed  to  be  on  the  order  of  the  thickness  of  the  plate,  while  the  tangential 
displacements  are  assumed  infinitesimal.  Finally,  KirchofP s  hypothesis  is  employed  with  tractions  on 
surfaces  paralle  to  the  middle  surface  assumed  negligible  and  strains  varying  linearly  with  the  plate 
thickness. 

Using  these  assumptions  the  governing  equations  for  the  plate  motion  may  be  written  as 


DV  w+m-r^=p,  +  n^  —  +  n^-z^  +  2n 


dx  dy 


dP 

=  0 


dx^ 


dr 


-^  +  ^  =  0 
dx  ^ 

where  the  two-dimensional  Laplacian  operator  and  the  internal  forces  are  defined  as 


n=B 


du  \(  awY 
rrz  +  -\  -zpr  1  +V\ 

cbc  2^aj: 


av  1  aw 
^^2  ¥ 


dv  1 

aj)  2 


+t) 


ax 


X  2(^ax 


(5-30-a) 

(5-30-b) 

(5-30-c) 

(5-31) 

(5-32-a) 

(5-32-b) 


B  l-u 


n^  =  - 


du  dv  dw  aw 
a^  dx  dx  dy 


(5-32-c) 


The  rigidity  constants  for  the  bending  and  membrane,  D  and  51 ,  and  the  mass  per  unit  area  are 
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D  = 


E.h^ 


o  Eh  _  .  f 


,  ,  „  ,  ,  ........  (5-33) 

12(1 -u")  (1-u  ) 

In  equations  (5-30)  through  (5-33),  u,  v ,  and  w  are  the  displacements  in  the  x ,  ,  and  z  directions 

respectively;  h  and  are,  respectively,  the  thickness  and  mass  density  of  the  plate;  is  external 
transverse  pressure  acted  on  the  plate. 

The  dynamic  equations  of  equilibrium  defined  above  are  dimensional.  To  make  them  compatible 
to  the  governing  equations  for  flow,  they  are  non-dimensionalized  with  respect  to  the  length  of  the 

plate.  Suppose  the  plate  is  square  and  the  length  is  / ,  the  non-dimensional  parameters  are  defined  as 

_  Pz 


A  y 

u  V  w  ,  n 

u  =  -,  v  =  -,  w  =  —  ,  n  =  -j,  Pz  -  ,  - 2 
till 


(5-34) 


Using  equation  (5-34),  the  non-dimensional  form  of  the  dynamic  equations  of  the  plate  becomes 

(5-35-a) 


dn^ 

O 

li 

A 

ox 

dy 

n„ 

xy  1 

-^  =  0 

9x 

dy 

and  the  two-dimensional  Laplacian  operator  and  the  internal  forces  are  defined  as 


=  t„l 


9m  if  dw\ 

— -1--  —  +v 

dx  2  9x  J 


ny=tj 


9v  1 

dy  2 


9y 

V  y 


9v  1 
9y  2 


9^ 


v^3^yj 


9m  1[^9w^ 

[^9x  2^^  9r  J  J 


1-u  f 9m  ,  9v  _  9w9w'^ 
“  2  ac  9y  j 


(5-35-b) 

(5-35-c) 

(5-36) 

(5-37-a) 

(5-37-b) 

(5-37-c) 


The  nondimensional  dynamic  pressure  X ,  mass  ratio  fJ,^ ,  and  internal  force  coefficient  are 
D  p^h  X[h) 


132 


(5-38) 


Computational  Mechanics  Laboratory,  University  of  Arkansas  at  Fayetteville  Selvam,  Qu  and  Zheng 


5.5  Numerical  Procedure 


5.5.1  Aerodynamic 

Solutions  to  equation  (5-10)  are  obtained  numerically  using  the  implicit  approximately  factored 
finite  difference  algorithm  of  Beam  and  Warmin^^,  employing  a  Newton  like  subiteration 
procedure^’*.  The  numerical  algorithm  is  obtained  from  equation  (5-10)  by  utilizing  either  a  two-  or 
three-point  backward  time  differencing  and  linearizing  about  the  solution  at  subiteration  level  .  The 
choice  of  the  first  or  second-order  temporal  accuracy  is  retained  in  the  following  iterative  ^proach  by 
specifying  ^  =  or  ^  =  1/2 ,  respectively.  The  nxunerical  algorithm  is  written  in  an  approximately 
factored,  delta  form  as 


(j  +<l>‘At^5^ 
(j-'f'  +(l>‘At,5^ 


-U" 


LV^/ 


!Ll 


dF” 

1 

dPf 

dU 

Re, 

du 

dH” 

1 

dH^ 

du 

Re, 

du 

^  .1 

-1/7+1 

yj 


1  dd;  ^ 


du  Re,  du 


jp+i 


MJ 


.  1  . 

Re^ 


+  5„ 


O’ --^a; 

Re, 


h'’-4-h: 

Re, 

(5-39)^ 


where 


f  =  — ,  AU  =  U'’^'-U‘’ 

l  +  (p 

Here,  U’’  is  the  subiteration  approximation  of  .  Therefore,  U’’  =  U‘  for 


(540) 


=  and 


UP  as  ^  oo . 

It  should  be  noted  that  with  this  subiteration  approach  the  right-hand  side  of  equation  (5-39) 
represents  the  numerical  approximation  of  the  governing  equation,  while  the  left-hand  side  vanishes  as 
—>  oo .  The  left-hand  side,  therefore,  may  be  modified  without  loss  of  formal  accuracy  provided  a 

sufficient  number  of  subiterations  are  employed.  In  particular,  a  time  step  on  the  left-hand  side  of  the 
equation,  At^ ,  may  be  chosen  independently  from  the  physical  time  step  At  on  the  right-hand  side  to 
enhance  the  stability.  Also,  the  right-hand  side  of  equation  (5-39)  can  be  modified  to  include  a  higher 
order  upwind  algorithm,  lagged  boundary  conditions  or  lagged  k  —  e  turbulence  modeling  without 


133 


CHAPTER  5:  Three  Dimensional  Aeroelastic  Solver  for  Nonlinear  Panel  Flutter 


destroying  the  implicit  nature  of  the  algorithm.  Left-hand  side  efficiency  improvements  can  also  be 
implemented.  The  numerical  procedure  has  been  modified  to  include  diagonalization,  following  the 
approach  by  Pulliam  and  Chaussee^®.  Although  the  diagonalized  form  of  the  ADI  scheme  is  only  first- 
order  time-accurate,  when  coupled  with  subiterations,  higher  order  time  accuracy  may  be  retained.  The 
numerical  scheme  (5-39)  reverts  to  the  standard  first  order  Beam-Warming  procedure  for  ^  =  , 

At^  =  At  ,aad  =  . 

In  equation  (5-39)  all  spatial  derivatives  are  approximated  by  the  second-order  accurate 
differences,  and  common  forms  of  both  implicit  and  explicit  nonlinear  dissipation^’  are  employed  in 
order  to  preserve  numerical  stability.  The  gird  speeds  x, ,  ,  and  z,  are  computed  in  a  manner 

consistent  with  the  temporal  derivatives  of  the  conserved  variables  in  equation  (5-39). 


5.5.2  Structure 

The  structural  equations  (5-30)  are  discretized  using  a  finite  difference  procediure  in  space  and  the 
Newmark-)3  method  in  time. 

The  Newmark-  integration  method  is  applied  to  solve  the  dynamic  equilibrium  of  motion.  It  is  a 
second-order  accurate  in  time  and  unconditionally  stable  approach.  In  this  method  the  following 
assumptions  for  the  velocity  and  displacement  in  the  time  interval  to  f  +  Af  are  employed: 

-1-  (1  -  d)w,  +  At  (5-41a) 

+  0/2  -  Ot)w,  +  aw)+^  At  ^  (5-41b) 

where  a  and  P  are  integration  parameters  that  determine  the  stability  and  accuracy  of  the  method. 
For  an  imconditionally  stable  integration  scheme,  Newmark  originally  proposed  that  a  =  1/4  and 
5  =1/2,  in  which  case  equation  (5-41)  corresponds  to  the  constant  average  acceleration  method. 

From  equation  (5-41),  may  be  solved  in  terms  of  as 

=  «o  -a^w,-a,w,  (5-42a) 

Substituting  equation  (5-42a)  for  in  equation  (5-41a)  yields 


+ 


or  in  a  simple  form 

<A,=«i  -a^w,-a^w, 


(5-42b) 
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The  constants  in  equation  (5-42)  are  defined  as 


_  1  _  8  _  1  _r  1  1 

a{Aty  ’  oAt  ’  oAt  ’  .2a 


_8  .  _At(5  . 

^4  —  f  5  ^5  —  ^  ^ 

a  2  a 


(5-43) 


Both  the  acceleration  w,^.^  and  velocity  at  time  t  + At  are  expressed  in  terms  of  the  unknown 
displacement  at  time  t  + At  and  the  parameters  in  time  .  Using  equation  (5-42),  the  governing 
equation  (5-35a)  of  the  plate  can  be  discretized  in  time  domain  as 

+  pI^  (5^) 

where  the  equivalent  external  force  is  defined  as 


PnA,  =  +  "2^/  + 


(5-45) 


Equations  (5-35b)  and  (5-35c)  are  static  equilibrium,  there  is  no  any  change  for  them. 

The  standard  five-  and  three-point  centered  finite  difference  expressions^^  are  used  to  approximate 
the  foiufh-,  second,  and  first-order  spatial  derivatives  in  equations  (5-44),  (5-35b),  (5-35c),  and  (5-37). 
The  derivatives  expressions  are  given  by 

j  «  -4m'.^,,,  +6w,.^.  -4w,_,,^.  + 

1  /  .  c  /I  _i_  ) 

^  1  r  /  \ 

I  ^  {AxA^  ~ 

^  d^w  ^  1  /  \ 

^  .*44x4;™"' 


3w^  1 


dx  \j  2Ax 


(5-46) 
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For  the  fourth-order  expressions,  i  =  +  2  to  ipe~^  (where  ip^  and  z^^  denote  the  beginning  and 

end  of  the  plate  in  the  x  direction)  and  j  =  jp^  +  2  to  jp^  -  2  (where  jp^  and  jp^  denote  the 

beginning  and  end  of  the  plate  in  the  direction).  For  z  =  z^^  + 1  or  z^^  —  1  and  j  —  j p^  + 1  or  y  1 , 

these  equations  are  modified  by  inserting  the  appropriate  boundary  condition  information. 

Introducing  the  derivative  expressions  in  equation  (5-46)  into  the  equivalent  static  equilibrium  (5- 
44),  the  finite  difference  expression  of  the  dynamic  equation  (5-35a)  for  AY  =  4y=  without 


boundary  condition  is  obtained  as 


(5-47a) 

where  is  replaced  by  jU  for  simplicity.  Similarly,  the  finite  expressions  for  the  govenung 

equations  (5-35b)  and  (5-35c)  are  given  by 

C„,  u{i-\J)  +  u{i+\,j)  +C„2  «(z,y-l)  +  tz(z,y  +  l)  +A„zz(z,y) 

+  C,[v(z  - 1,  y  - 1)  -  v(z  + 1,;  - 1) + v(z  +  l,y  + 1)  -  v(z  - 1,  y  + 1)] 

+  Qi [Mi  - 1. 7  - 1)  -  Mi  + 7  - 1)  +  Mi  +  l.y  + 1)  -  Mi  -hj+ 1)]  (547b) 

[->v(z,y-l)  +  M<z,y>l)] 

+  C„2  [-  Mi  ~  1>7)  +  Mi  + j)\Mi  - J)  -  j)  +  Mi  + )] 

+  Mi  - 1.  j)  +  Mi  +  hjilMij'  - 1)  -  '2'Mij')  +  Mi^  J  + 1)1= ^ 
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C„2  vii-lj)  +  v(i+lj)  +C„,  v(/J-l)  +  v(i,7  +  l)  +AXiJ) 

+  C^[u(i  - 1, 7  - 1) -  u(i + 1,7  - 1) + m(/ ■+ 1,7  + 1) -  u{i  - 1, 7  + 1)] 

+  C^,  [m</  - 1, 7  - 1)  -  w(i  + 1, 7  - 1)  +  w(i + 1,7  + 1)  -  w(/  - 1,7  + 1)] 
[-w(i-l,7)  +  M<i  +  l,7)] 

+  C„2  [-  w(i,  7  - 1)  +  Mi,  7  +  7  - 1)  -  ^Mi,  j)  +  MU  + 1)] 

+  Mi,  7  “  1)  +  Mi,  7  +  -hj)-  '^Mi,j)  +  Mi  + 1. 7')]  =  0 


(5-47c) 


where. 


r  -  h  r  r  r  c  =-^  C 

2h^  ’  8/j^  ’  16A^  ’  2h^'  4h^ 

(5-48a) 

4=-2C„,  +  Q3  (5-48b) 

After  the  proper  boundary  conditions  are  applied  to  Equation  (5-47),  they  can  be  solved.  Assembling 
equation  (5-47),  the  governing  equations  of  the  discrete  model  of  the  large  deflection  plate  at  the 
current  time  step  is  obtained  as 

A{u,  V,  w)  =  B(u,  V,  >v)„  (5-49a) 

A(u,v,w)^  w  =  B(u,v,  w)„  (5-4%) 

/1(m,v,w)„m' =  5(m,v,>v)^  (5-49c) 

Obviously,  the  coefficient  matrices  and  vectors  on  the  right-hand  sides  of  equation  (5-49)  are  the 
functions  of  the  unknowns.  This  means  equation  (5-49)  is  nonlinear.  Iterations,  not  only  within 
equations  of  the  equation  (5-49)  but  also  among  them,  are  necessary  to  solve  for  the  unknowns  from 
this  equation. 


5.5.3  Fluid-Structure  Interaction 

The  fluid-structure  interaction  analysis  may  be  performed  using  the  following  three  main  steps  in 
each  time  step.  (1)  Calculate  the  responses  of  the  panel  under  the  fluid  pressure  and  other  external 
forces;  (2)  Move  the  fluid  grid  using  the  displacements  from  the  panel  analysis;  (3)  Perform  the  fluid 
analysis  imder  the  moved  grid  and  compute  the  pressure  on  the  surface  of  the  panel  due  to  the  its 
movement.  To  start  the  numerical  sunulation,  some  initial  perturbation,  initial  displacements,  velocity, 
or  accelerations,  is  usually  assigned  to  the  panel.  Steps  one  and  three  are  corresponding  to  structural 
and  fluid  analysis  respectively.  For  each  step,  iterations  are  required  because  the  governing  equations 
of  the  structure  and  fluid  are  nonlinear. 
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The  pressure  resulted  from  the  fluid  side  is  required  when  calculating  the  responses  of  the  panel. 
On  the  other  hand,  the  displacements  of  the  panel  should  be  available  before  performing  the  analysis 
of  the  fluid  and  obtaining  the  pressure  to  be  assigned  on  the  panel.  Clearly,  the  fluid  and  the  panel  are 
coupled.  If  only  one  cycle,  from  sfructure  to  fluid,  is  used  at  each  time  step,  big  errors  will  sometimes 
be  resulted  due  to  the  lagged  fluid-stracture  coupling.  This  error  will  lead  numerical  instability  when  a 
relatively  long  time  simulation  is  performed.  To  eliminate  the  possible  numerical  instability,  the 
subiteration  between  the  fluid  and  stracture  is  implemented. 

As  we  know,  the  computational  effort  of  the  fluid  is  much  more  expensive  than  the  stracture,  the 
following  main  logic  is  used  to  reduce  the  computational  work. 

1.  Assign  the  initial  state  of  the  panel; 

2.  For  each  time  step,  do: 

2.1.  For  each  subiteration  p,  do: 

2.1.1  Evaluate  the  equivalent  external  forces  acting  on  the  plate  using 
P't+^  =  ^p’’  +  a^w,  +  w,  +  OjVt;  and  solve  equation  (5-49)  iteratively; 

2.1.2  Move  the  gird  using  the  displacement  vector  «,  v,and  w  of  the  panel; 

2.1.3  Perform  the  fluid  analysis  and  find  the  pressure  acting  on  the  panel; 

2.1.4  Set  =  +  and  go  back  to  step  2.1.1; 

2.2.  Evaluate  the  velocities  and  accelerations  at  current  time  using  equation  (5-42); 

2.3.  Ouq)ut  the  results; 

3.  Output  the  results. 

5.6  Results  and  Discussions 

The  problem  to  be  investigated  is  the  flow,  inviscid  flow  at  this  time,  over  a  two-dimensional 
square  flexible  panel  of  length  /,  as  shown  in  Figure  5.1.  Unless  otherwise  noted  the  panel  has  the 
following  properties:  thickness  h/l  =  0.002 ,  mass  ratio  p=  and  Poisson’s  ratio  V  =  0.3 .  For  all 
cases  fieestream  pressure,  ,  is  specified  on  the  underside  of  the  panel.  The  coordinate  of  the  flow 
-  .Yvz  and  are  shown  in  Figure  5.1. 

Numerical  boundary  conditions  for  the  panel  provide  the  connection  between  the  aerodynamic  and 
structural  equations.  For  viscous  flow  computations,  a  no  slip  condition  is  applied  on  the  plate  surface. 
This  requires 
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=  '’y=yb 

where  Xf,  and  denote  the  velocity  of  the  moving  boundary  with  =;)>*=  0  in  the  static  case.  The 
remaining  two  conditions  are  the  adiabatic  wall  condition  and  the  normal  momentum  equation: 


dri 


Along  the  inflow  boundary,  all  dependent  variables  are  assigned  their  respective  freestream  values  for 
supersonic  flow,  whereas  for  the  subsonic  flow  characteristic  boundary  conditions^*  are  applied.  On 
the  top  boundary,  either  extrapolation  or  characteristic  conditions  are  specified  for  supersomc  or 
subsonic  flows  respectively.  On  the  outflow  boundary,  first-order  accurate  extrapolation  of  the 
dependent  variables  is  employed  in  all  cases,  corresponding  to  the  condition 


For  the  inviscid  case,  the  boundary  conditions  along  the  plate  are  modified  by  setting  the  fluid 
velocity  component  normal  to  the  surface  equal  to  the  corresponding  values  for  the  plate.  Finally,  a 
slip  condition  is  implemented  by  using  the  second-order  extrapolation  for  the  tangential  velocity 


component. 


Figure  5.1.  Geometry  of  the  flow  and  panel 
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Boundary  conditions  for  the  plate  are  specified  for  simply  supported  edges.  No  deflection  is 
allowed  along  the  edges  of  the  plate,  i.e.,  m  =  v  =  w  =  0.  No  moment  is  on  the  edges,  that  is, 
d^wfdn^  =  0. 

The  grid  used  for  the  panel  is  20  X  20  with  equal  space.  The  finite  difference  grid  for  the  flow  is 

xXyXz  =  4  x4  x2  as  shown  in  Figure  5.2.  The  minimum  space  in  the  z-direction  is  0.001. 


5.6.1 .  Amplitudes  of  Limit  Cycle  Oscillation  (LCO) 
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As  mentioned  above,  panel  flutter  is  induced  by  supersonic  air  flow  on  a  panel.  The  aerodynarmc 
pressure  acting  on  the  panel  is  a  function  of  the  panel  motion  itself.  Therefore,  it  is  very  essential  to 
accurately  predict  the  pressure.  Even  though  there  are  several  aerodynamic  theories,  most  of  them  are 

V 

linearized,  approximate.  For  the  current  research,  the  aerodynamic  loads  are  resulted  from  the 
simulation  of  the  Euler  flow  on  the  surface  of  the  nonlinear  panel.  The  amplitudes  of  the  limit  cycle 
oscillation  resulted  from  the  present  computation  will  be  compared  with  those  from  approximate 
aerodynamic  theories.  Two  of  them,  folly  linearized  potential  flow  and  qusi-steady  piston  theory,  will 
considered  at  this  time. 

The  folly  linearized  aerodynamic  theory  was  used  to  estimate  the  aerodynamic  pressure  on  the 
panel  by  Cunningham^*  and  Dowell*^  about  two  decades  ago.  This  theory  is  based  on  the  inviscid, 
potential  flow.  It  is  usually  valid  for  a  low  Mach  number  (AC  =!)•  The  amplitudes  of  the  LCO  at 
(0.75,0.50)  in  flie  panel  coordinate  obtained  from  this  theory  and  the  present  calculation  are  plotted  in 
Figure  5.3.  Only  three  Mach  numbers,  1.2,  1.414,  and  1.6,  are  compared  due  to  the  availability  of  the 
former  results. 

As  shown  in  Figure  5.3,  the  critical  flutter  dynamic  pressures,  A ,  are  216,  462,  608  for  the  Euler 
flow  while  they  are  191, 421,  and  591  for  the  linearized  potential  flow.  The  former  pressures  are  13%, 
9%,  and  3%  higher  than  the  later  respectively.  The  present  amplitudes  are  lower  than  those  from  the 
linearized  potential  flow  theory  as  shown  in  Figure  5.3.  With  the  increase  of  the  Mach  number,  the 
difference  becomes  smaller  and  smaller.  The  reason  might  be  that  die  nonlinearity  is  ignored  in  the 
linearized  potential  theory.  The  percent  difference  of  the  amplitudes  resulted  from  the  potential  flow 
with  respect  to  the  results  from  the  Euler  flow  are  listed  in  Table  5.1.  Generally,  the  difference  at 
higher  Mach  number  is  much  smaller  than  the  lower  Mach  number. 
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Figure  5.3.  Comparison  of  LCO  amplitudes:  Potential  flow  and  Euler  flow 


Table  5.1.  Percent  differences  of  the  LCO  amplitudes 


A 

Diff. 

A 

Diff. 

1.2 

250 

75.3 

1.6 

650 

24.9 

1.2 

300 

62.0 

1.6 

700 

14.1 

1.414 

500 

88.7 

1.6 

750 

13.0 

1.414 

550 

60.8 

1.6 

800 

16.5 

1.414 

600 

1.6 

18.3 

The  qusi-steady  first  order  piston  theory  was  proposed  to  predict  the  aerodynamic  pressure  by 
Ashley  and  Zartarian  in  1956.  Since  the  three  dimensionality  and  the  unsteadiness  of  the  air  flow  was 
ignored  in  this  theory,  it  cannot  be  applied  for  the  air  flow  with  Mach  number  close  to  one.  It  mostly 

applied  for  the  airflow  with  large  Mach  number  (M^>  ^/2  ).  The  LCO  amplitudes  obtained  from  the 
two  theories  are  shown  in  Figure  5.4.  The  corresponding  difference  is  Usted  in  Table  5.2.  It  is  shown 
clearly  in  Figure  5.4  that  the  present  amplitudes  are  higher  than  those  from  the  qusi-steady  piston 
theory.  The  percent  difference  reduces  with  the  increase  of  the  dimensionless  pressure.  The  difference 
for  Mach  number  1.6,  for  example,  reduces  from  44.4%  at  A  =  650  to  9.3%  at  A  =  850 .  This  also 
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happens  for  Mach  number  2.0.  The  most  interesting  phenomenon  is  that  the  shapes  of  curve  resulted 
from  the  two  theories  are  almost  same  for  Mach  number  1.6  and  2.0  respectively.  Therefore,  there  is 
only  some  shifting  difference. 


Figure  5.4.  Comparison  of  LCO  amplitudes:  Piston  theory  and  Euler  flow 


Table  5.2.  Percent  differences  of  the  LCO  amplitudes 


A 

Diff. 

X 

Diff. 

1.6 

650 

44.4 

2.0 

900 

41.8 

1.6 

700 

19.8 

2.0 

950 

17.0 

1.6 

750 

14.0 

2.0 

1000 

10.0 

1.6 

800 

11.2 

2.0 

1050 

7.4 

1.6 

850 

9.3 

2.0 

1100 

6.2 

As  mentioned  in  the  introduction,  Gordnier  and  Visbal  have  also  used  Euler  equations  to  predict 
the  aerodynamic  loads  very  recently.  The  present  amplitudes  are  compared  with  theirs  in  Figure  5.5. 
Only  the  results  at  Mach  number  1.2  are  shovra  due  to  the  availability  of  the  later  results.  They  are 
very  close  as  expected.  However,  there  is  a  big  difference  for  the  critical  dynanic  pressure.  The 
prediction  of  the  critical  dynamic  pressure  is  very  computationally  expensive  when  the  high  fidelity 
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equations  are  used  to  compute  the  aerodynamic  loads.  Since  the  critical  pressure  can  not  be  predict 
explicitly,  six  to  ten  simulations  for  the  dimensionless  pressures  around  the  critical  are  usually 
required  even  though  the  linear  interpolation  searching  technique  is  employed.  For  example,  about  one 
week  was  used  to  find  the  critical  pressure  when  the  error  is  required  to  be  less  than  1.  The 
computation  was  performed  on  a  Sun  4500  computer  with  four  of  its  eight  processors. 


Figure  5.5.  Comparison  of  LCO  amplitudes:  Euler  flow 

The  deflections  of  the  panel  at  the  end  of  the  500  dimensionless  time  are  plotted  in  Figures  5.6  and 
5.7.  For  these  figures,  Mach  number  M„=1.6,  dimensionless  dynamic  pressure  A  =  650  and 
A  =900  are  used.  The  following  three  phenomina  may  be  easily  obtained  firom  these  plots.  (1)  The 
maximum  deflection  is  approximately  located  on  the  right  one  quarter,  i.e.,  x^jl  =0.75  and 

j;^//  =0.5  .  This  is  shown  clearly  in  Figures  5.6(a)  and  5.7(a).  (2)  The  deformations  are  symmetric 
with  respect  the  central  line  of  the  panel  in  the  x  direction,  that  is,  =  0.5 .  (3)  In  the  y  direction, 
the  shape  of  the  deformation  is  very  close  to  the  first  modeshpe  of  the  panel  itself  These  three 
phenomena  are  very  similar  to  those  based  on  the  approximate  theories. 
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Figure  5.8.  Surface  pressure  contour  of  the  panel  for  ==  1 .6 :  (a)  A  =  650 ,  (b)  A  =  900 
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The  pressures  and  the  flow  velocity  in  the  x  -direction  on  the  surface  of  the  panel  at  the  same  time 
mentioned  above  are  shown  in  Figures  5.8  and  5.9  respectively.  Clearly,  they  are  symmetrical  with 
respect  to  the  central  line  of  the  panel  in  the  x -direction,  that  is  =  5  .  As  we  know,  the 

dimenionless  air  pressiure  of  the  freestream  is  0.2790  as  shown  in  the  left  one  fifth  part  in  Figure  5.8. 
Usually,  there  are  three  net  pressure  ranges.  They  approximately  locate  on  the  left  one  quarter,  middle 
one  half,  and  right  one  quarter  of  the  panel  respectively  as  shown  in  Figure  5.8.  For  convenience,  they 
are  called  range  I,  II,  and  HI  from  lest  to  right.  Since  the  net  pressures  in  range  I  are  very  close  to  the 
reference  pressiue  0.2790,  they  are  not  clearly  shown  in  Figure  5.8.  The  net  pressure  in  range  HI  is 
much  higher  than  that  in  range  II  which  is  much  higher  than  that  in  range  I.  This  is  the  reason  that  the 
maYimiim  deflection  of  the  panel  is  always  very  close  the  three-quarter  of  the  panel  from  left. 
Furthermore,  the  pressure  within  ranges  II  and  ni  usually  have  different  sign.  Of  course,  there  is  a 
pressure  shock  between  these  two  ranges.  The  phenomenon  of  the  distribution  of  the  pressure  may  be 
explained  from  the  air  velocity  distribution  on  the  surface  of  the  panel  shown  in  Figure  5.9. 


5.6.2  Flutter  Frequency 

The  LCO  in  time  domain  is  plotted  for  =1.2  and  A  =300  in  Figure  5.10(a).  The  FFT 
transformation  of  the  oscillation  in  frequency  domain  is  shown  in  Figure  5.10(b).  Initially,  the 
amplitude  grows  rapidly  due  to  the  negative  damping  from  the  aerodynamic  pressure.  With  the 
increase  of  the  amplitude  the  membrane  tensile  forces  become  significant  which  bound  the  flutter 
amplitude  at  a  constant  some  time  later.  It  can  be  foimd  from  Fig.  5.10(b)  that  the  fundamental  flutter 
frequency  of  the  LCO  is  0.0854. 

As  we  know,  the  natural  circular  frequencies  of  the  linear  simply  supported  thin  plate  can  be 
expressed  as^"* 

where  m  =  •  •  • .  The  corresponding  nondimensional  frequencies  are 


/>  ^  (  2  2\  I  ^ 

J  mn  /^rr  7  '  M  ^  7 


2VJ 


Psh 


(5-51) 


Using  the  dimensionless  dynamic  pressure  A  and  mass  ratio  ,  equation  (5-51)  may  be  rewritten  as 


fm 


(5-52) 
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Similarly,  the  flutter  frequencies  and  the  corresponding  natural  frequencies  may  be  calculated. 
They  are  plotted  in  Fig.  5.11.  In  these  pictures,  “Linearl”  and  “Linear2”  denote  the  first  and  second 
natural  frequencies  respectively.  The  natural  frequencies  decrease  with  the  increase  of  the 
nondimensional  dynamic  pressure  which  can  be  explained  clearly  from  equation  (5-52).  As  we  know, 
the  panel  becomes  more  flexible  when  the  dynamic  pressure  increases.  Under  the  same  air  flow,  the 
responses  should  be  bigger  for  the  higher  dynamic  pressure.  The  nonlinear  effect,  therefore,  gets  more 
significant.  Consequently,  the  flutter  frequencies  obtained  from  both  Euler  flow  theory  and  piston 
theory  increase  very  slightly  when  the  dynamic  pressure  increases  even  though  the  corresponding 
natural  frequencies  decrease.  For  the  higher  Mach  numbers,  M„  =1.6  and  =  2.0 ,  the  flutter 
frequencies  are  almost  constant.  The  flutter  frequencies  resulted  from  the  present  simulation  are  a  little 
lower  than  those  from  piston  theory.  Because  the  flutter  frequency  is  pertaining  to  the  nonlinear 
model,  it  should  be  bigger  than  the  corresponding  natural  frequency.  From  Fig.  5.11  we  know,  almost 
all  the  flutter  frequencies  are  located  between  the  first  and  the  second  natural  frequencies  of  the 
corresponding  linear  models.  This  means  that  the  lowest  two  modes,  especially  the  first  mode,  have 
significant  contribution  to  the  LCO. 


(a) 
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Figure  5.1 1.  Comparison  of  the  frequencies: 

(a)  =  1.2;  (b)  =  1.414;  (c)  =1.6;  (d)  =2.0 


5.6.3.  Effect  of  the  Panel  Thickness 

In  the  following  discussion,  we  assume  the  nondimensional  dynamic  pressure  A ,  mass  ratio  , 

and  the  parameters  of  the  air  flow  are  constants  when  the  thickness  h  changes.  This  means  that  the 
panel  and  the  air  flow  do  not  change  in  the  nondimensional  sense  for  different  thickness. 


5.6.3.1  Linear  Case 

The  linear  governing  equations  of  the  panel  under  the  Euler  flow  can  be  easily  obtained  from 
equation  (5-35a)  as 


u  4  d^w 


(5-53) 
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(a) 


(b) 

Figure  5.12.  Linear  deflection  of  the  panel:  (a)  A  =  650 ;  (b)  A  =  900 
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As  we  know,  the  aerodynamic  pressure  resulted  from  the  Euler  equations  is  only  dependent  on  the 
deflection  w  of  the  panel  during  the  simulation.  Hence,  it  will  not  be  affected  by  the  thickness.  The 
items  on  the  left-hand  side  also  have  nothing  to  do  with  the  thickness.  Therefore,  the  deflection  of  the 
panel  is  supposed  to  be  independent  of  the  thickness.  Since  the  critical  dimensionless  dynamic 
pressure,  ,  can  be  determined  from  the  linear  panel,  it  is  also  independent  of  the  thickness. 

The  deflections  of  the  panel  for  =  1.6 ,  X  =  650  and  A  =  900  are  plotted  in  Figure  5.12.  In 

these  results,  three  different  thicknesses,  0.02,  0.002,  and  0.0002,  are  used.  The  deformations  blow  up 
at  time  90.58  in  Figure  5.12(b)  due  to  the  unreasonable  deformation  of  the  panel.  The  deformations  are 
exactly  same  for  different  thicknesses  of  the  panel. 

5.6.3.2  Nonlinear  Case 

For  the  nonlinear  case,  the  effect  of  the  thickness  becomes  a  little  complex.  Let’s  discuss  the 
influence  of  the  thickness  on  the  deformation  w  at  first.  Introducing  equation  (5-37)  into  equations  (5- 
35b)  and  (5-35c)  leads 

d^u  f  d^v  dw  d^w  ^  / d^u  d\  9w  d^w  ^  ^  ^ 

dx^^dxdx^  ^ ^dxdy  dy  dxdy  ^  dy  dxdy  ^ 

(5-54-a) 

d^u  /,  /  d^u  ^  3^v  ^  3w  3^w  ^  3w3^w 

dy^  ^  dy^  3x3y  dx  dxdy^  |^3x3y  3x^  dx  dxdy  3>’  dx^  ^ 

(5-54-b) 

These  two  equations  are  independent  of  the  thickness.  Therefore,  only  Equation  (5-35a)  depends  on 
the  thickness  through  the  internal  forces  expressed  in  equation  (5-37).  As  shown  in  equation  (5-37), 
the  internal  forces  reduce  with  the  increase  of  the  thickness.  This  leads  the  decrease  of  the  nonlinear 
effect.  Hence,  the  deflection  w  will  increase.  This  is  shown  clearly  in  Figures  5.13  and  5.14.  The 
deflection  of  the  case  h  =  0.02  is  much  bigger  than  the  case  h  =  0.002  and  the  later  is  much  bigger 
than  the  case  h  =  0.0002  .As  shown  in  Figures  5.13(b)  and  5.14(b),  the  deflections  are  very  close  to 
each  other  for  different  thicknesses  at  beginning.  The  reason  is  that  these  deflections  are  very  small 
and  the  nonlinear  effect  is  unimportant.  They  are  close  to  the  linear  case. 
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(a) 


(b) 


Figure  5.14.  Nonlinear  deflection  of  the  panel  for  A  =  900 
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Defining  the  ratio  of  the  deflection  w  with  respect  to  the  thickness  h  as  w-  w/h .  Substituting 
w  for  w  in  equations  (5-35a)  leads 


X  '  '  h 

in  which,  the  internal  forces  are 


12 


12 


diu  Bw'f 


2 


ydx  j 


3v 


dy"^  2 


+1) 


+U 


^  +  i 

dy  2 


ay 


dxdy 


(5-55) 


du  i^awY 
dx  2(^aAi:  J 


6  1-v 


■  X 


du  ^  dv  ^  awdyv^ 


(5-56-a) 


(5-56-b) 


(5-56-c) 


^  dx  dx  dy  J 

where,  u  =  ujh^  and  v  =  v//?^ .  Using  equations  (5-56),  the  governing  equations  (5-35b)  and  (5-35c) 
become 


d^u  dw  d^w 
dx^  dx  dx^ 


(  :j2 


d^v  dw  d^w  ^ 


dxdy  dy  dx^ 


-(l-u] 


di^u  d^v  dwd^w  dw  d^w 


dy^  ^  axa>'  ^  dx  dy^  ^  a>'  dxdy 

< 

(5-57-a) 


=  0 


a^v  dw  a^w 

ay  dy  a^^ 


/  :»2 


-i-u 


a^M  aiv  a^w 

+■ 


ardj  dx  dxdy 


+  (1 


/a^M 

■^|axaj 


a^v  dw  a^w  awd^w^ 

dxajj'  ^  dx^  ^  dx  dxdy  dy  dx^ 


:0 


(5-57-b) 


The  ratios  of  the  deflection  for  A  =  650  and  A  =  900  are  shown  in  Figures  5.15  and  5.16 
respectively.  It  can  be  seen  that  the  amplitudes  of  the  ratio  become  stable  after  some  time,  200  for 
example.  The  stable  amplitudes  are  very  close  for  the  same  nondimensional  dynamic  pressure.  This 
means  that  the  ratio  of  the  LCD  amplitude  is  independent  of  the  thickness.  We  also  can  conclude  that 
the  ratio  of  pjh  has  nothing  to  do  with  the  thickness  for  the  Mach  number. 
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This  phenomenon  can  be  explained  clearly  if  the  qusi-steady  piston  aerodynamic  theory  is 
considered.  In  this  theory,  the  aerodynamic  pressure  acting  on  the  panel  is  assumed  as 


r  z 


Ml  -1  dt  ^ 


Its  nondimensional  form  is  given  by 

/ 


■JmI-1 


dx^  Ml -\  dt 


V 


Introducing  equation  (5-59)  into  equation  (5-55)  results 


a'w 


MsP. 

^dw  Ml—2dw^ 
_ 1 _ 22 - 

^Ml-l 

dx  Ml-l  dt 

K  ^  y 

a^w _  d^w 


dx 


w 

T  +  Ms»y 


dy 


■  +  2ju,n^ 


dxdy 


(5-58) 


(5-59) 


Clearly,  the  above  equation  has  nothing  to  do  with  the  thickness.  Actually,  this  is  true  for  supersonic 
flow. 

Although  the  thickness  does  not  affect  the  ratio  of  the  LCO  amplitude  after  the  oscillation  become 
stable,  it  does  affect  the  way  to  the  LCO.  For  example,  the  ratio  of  the  amplitude  for  h  =  0.02 
increases  very  slowly  as  shown  in  Figure  5.15(a).  It  takes  about  1 10  nondimensional  time  for  the  ratio 
to  close  to  the  stable  value.  When  the  thickness  decreases  by  one  tenth  and  one  hundredth,  the 
computer  times  reduce  to  about  80  and  10,  respectively,  as  shown  in  Figure  5.15.  As  we  know,  the 
deflection  is  very  small  at  beginning  and  fall  into  the  linear  or  weak  nonlinear  range.  According  to  the 
results  in  Figures  5.13  and  5.14,  the  amplitudes  should  be  very  close.  Therefore,  the  ratio  which  is 
scaled  by  the  inversion  of  the  thickness  decreases  with  the  increase  of  the  thickness.  At  the  very 
beginning,  the  difference  is  very  close  to  1/h.  With  the  increase  of  the  amplitude,  the  difference 
decreases  very  fast.  After  the  oscillation  becomes  stable,  the  difference  becomes  insignificant. 


5.7  Summary 

The  three-dimensional  nonlinear  panel  flutter  at  supersonic  flow  is  analyzed  using  the  high-fidelity 
flow  equations,  Euler  equations.  These  flow  equations  are  solved  using  the  Beam-Warming,  alternate- 
direction,  implicit  scheme.  The  governing  equations  of  the  large  deflection  panel  are  based  on  the  von- 
Karman  theory.  Finite  difference  method  is  used  to  discretize  the  governing  equations  in  space  and 
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Newmark-Beta  integration  scheme  is  applied  to  solve  them  in  time  domain.  The  Newton-like 
subiteration  is  implemented  to  eliminate  the  lagging  errors  associated  with  the  exchange  of  the 
pressure  and  deformations  between  the  fluid  and  panel  at  their  interface.  The  numerical  sunulation  is 
performed  on  a  simply  supported  square  panel.  The  results  at  the  supersonic  air  flow  are  compared 
with  those  from  potential  flow  theory  and  qusi-steady  supersonic  theory. 

(1)  The  present  LCO  amplitudes  are  smaller  than  those  resulted  from  the  linear  potential  flow 
theory.  Generally,  the  difference  at  higher  Mach  number  is  smaller  flian  the  lower  Mach  number. 

(2)  The  LCO  amplitudes  obtained  from  Euler  flow  are  bigger  than  those  obtained  from  qusi-steady 
piston  theory.  Similarly,  the  difference  reduces  with  the  increase  of  the  dynamic  pressure. 

(3)  The  current  results  are  very  close  to  those  from  Reference  14  except  for  the  critical 
nondimensional  dynamic  pressure  at  Mach  number  1.2. 

(4)  The  flutter  frequencies  increase  very  slightly  with  the  increase  of  the  dynamic  pressure.  For 
high  supersonic  low,  =  2.0  for  example,  these  frequencies  become  a  constant  wdthin  a  wide 
dynamic  pressure  range.  The  current  flutter  frequencies  are  a  little  lower  than  those  from  qusi-steady 
piston  theory.  The  first  mode  of  the  panel  has  the  significant  contribution  of  the  panel  flutter. 

(5)  Assume  the  nondimensional  dynamic  pressure  A ,  mass  ratio  ,  and  the  parameters  of  the  air 
flow  are  constants.  The  deflection  of  the  panel  for  the  linear  case,  and  hence  the  critical 
nondimensional  dynamic  pressure,  is  independent  of  the  dimensionless  thickness  of  the  panel.  For  the 
nonlinear  case  of  the  panel,  the  nondiemsional  thickness  does  not  affect  the  ratio  of  the  LCO 
amplitude.  However,  it  does  affect  the  way  to  the  LCO. 
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CONCLUSIONS 


A  three-dimensional  solver  for  the  nonlinear  panel  flutter  at  supersonic  flow  has  been  developed. 
The  solver  has  been  verified  using  several  approximate  aerodynamic  theories  such  as  the  potential 
flow  theory  and  the  qusi-steady  supersonic  theory  on  a  simply  supported  square  panel  at  the  supersonic 
flow.  The  complex  problems  like  wing  flutter  can  be  solved  with  the  development  of  the  current  work. 

Finite  element  method  is  the  most  commonly  used  approach  for  discretizing  the  nonlinear  panel. 
The  resulted  nonlinear  equations  of  equilibrium  may  be  solved  using  three  formulations.  Two  of  them, 
the  Total  Lagrangian  formulation  and  the  Co-Rotational  formulation,  are  efficient  for  the  nonlinear 
panel  and  have  been  discussed  in  details  in  Chapters  1  and  2  respectively.  Because  only  the  panels 
with  very  simple  geometry  are  considered  at  this  time,  the  finite  difference  method  rather  than  the 
finite  element  method  for  solving  the  governing  equations  of  the  nonlinear  panel  was  implemented 
into  the  solver  to  save  the  computer  time.  The  Newmark-  jS  integration  scheme  was  applied  to  solve 

them  in  time  domain. 

In  the  fluid  side,  the  high-fidelity  flow  equations  -  Euler  equations  was  used  in  this  solver.  The 
viscous  effect  has  not  been  included  at  this  time.  The  finite  volume  method  is  applied  to  discrete  the 
fluid  in  space.  A  three-dimensional  aeroelastic  solver  using  QUICK  scheme  has  been  developed  using 
the  procedure  in  Ch^ter  4.  Because  of  the  use  of  sequential  solver  for  a  highly  nonlinear  Euler  (NS) 
equations,  the  implicit  sequential  solver  becomes  overall  explicit  solver  even  after  subiteration.  Hence, 
a  Beam-Warming,  alternate-direction,  implicit  scheme  solver  has  been  developed  for  the  elastic  solver. 

An  explicit  scheme  for  the  interaction  of  the  fluid  and  stracture  was  used.  The  Newton-like 
subiteration  was  implemented  to  eliminate  the  lagging  errors  associated  with  the  exchange  of  the 
pressiire  and  deformations  between  the  fluid  and  panel  at  their  interface. 

The  present  amplitudes  of  the  limit  cycle  oscillation  of  a  simply  supported  square  panel  are 
smaller  than  those  resulted  fi-om  the  linear  potential  flow  and  greater  than  those  fi-om  qusirsteady 
piston  theory.  The  difference  becomes  smaller  when  the  Mach  number  increases.  The  flutter  frequency 
and  the  effects  of  the  panel  thickness  are  also  discussed  in  this  report. 

Further  work  of  this  solver  is  underway  and  listed  in  the  following: 

•  Solve  for  the  nonlinear  panel  using  the  finite  element  method  together  with  both  formulations; 

•  Include  the  viscous  effect  in  the  solver  (use  NS  equations); 
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•  Solve  for  the  panel  flutter  with  complex  conditions; 

•  Extended  present  solver  for  the  wing  flutter  problems. 
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